[Date Index]
[Thread Index]
[Author Index]
RE: Re: help in generating a gaussian random variable
*To*: mathgroup at smc.vnet.net
*Subject*: [mg35545] RE: [mg35502] Re: [mg35488] help in generating a gaussian random variable
*From*: "Wolf, Hartmut" <Hartmut.Wolf at t-systems.com>
*Date*: Thu, 18 Jul 2002 03:06:52 -0400 (EDT)
*Sender*: owner-wri-mathgroup at wolfram.com
> -----Original Message-----
> From: Wolf, Hartmut
> Sent: Wednesday, July 17, 2002 4:33 PM
> Cc: 'Salman Durrani'
> Subject: [mg35545] RE: [mg35502] Re: [mg35488] help in generating a gaussian
> random variable
>
>
>
> > -----Original Message-----
> > From: Wolf, Hartmut
> > Sent: Wednesday, July 17, 2002 12:46 PM
> > Cc: 'Salman Durrani'
> > Subject: [mg35545] RE: [mg35502] Re: [mg35488] help in generating a gaussian
> > random variable
> >
> >
> >
> > > -----Original Message-----
> > > From: BobHanlon at aol.com [mailto:BobHanlon at aol.com]
> > > Sent: Wednesday, July 17, 2002 8:09 AM
> > > Subject: [mg35545] [mg35502] Re: [mg35488] help in generating a
> > gaussian random
> > > variable
> > >
> > >
> > >
> > > In a message dated 7/16/02 5:50:43 AM,
> > dsalman at itee.uq.edu.au writes:
> > >
> > > >I need to generate a Gaussian random variable y having mean
> > > =0 and variance
> > > >=x.
> > > >
> > > >The variance x is itself a gaussian random variable having a
> > > known mean
> > > >and
> > > >variance.
> > > >e.g. mean of x =5;
> > > >variance of x = 10;
> > > >
> > > >Can anyone suggest how to use this information to generate y ?
> > > >
> > >
> > > Needs["Statistics`NormalDistribution`"];
> > >
> > > data=RandomArray[
> > > NormalDistribution[0, Random[
> > > NormalDistribution[5, Sqrt[ 10]]]], {100}];
> > >
> > >
> > > Bob Hanlon
> > > Chantilly, VA USA
> > >
> >
> > Bob,
> >
> > sorry, but I'm wary of your proposed solution. This is
> > because RandomArray is an optimising function to speed up
> > over Table Random.
> >
> > A look into the package Statistics`ContinuousDistributions`
> > reveals that the arguments to NormalDistribution are
> > evaluated only once. So the variance x is not a random
> > variable. A statistics test of
> >
> > data2 = Table[
> > Random[NormalDistribution[0,
> > Random[NormalDistribution[5, Sqrt[10]]]]], {100}];
> >
> > should reveal this. To make it simple for me: taking samples
> > of 10001 events respectively, and comparing the cumulants e.g.
> >
> > Block[{$DisplayFunction = Identity},
> > ListPlot[
> > Take[Transpose[{Sort[#1], Range[Length[#1]]/Length[#1]}],
> > {1, -1, 100}],
> > PlotJoined -> True,
> > PlotStyle -> #2]] & @@@ {{data, {}}, {data2,
> > Hue[.7]}} // Show
> >
> > shows a great discrepancy. Perhaps you like to draw multiple
> > samples of data (not necessarily of data2) to compare.
> >
> > This also shows that the desired distribution is not a normal
> > distribution (when x is a random variable).
> >
> >
> > Even more revealing is to compare the distributions by eye view
> >
> > GraphicsArray[
> > Block[{$DisplayFunction = Identity, step = 100, scale},
> > scale = step/Length[#1];
> > ListPlot[{Plus[##]/2, -scale/Subtract[##]} & @@@
> > Partition[Take[Sort[#1], {1, -1, step}], 2, 1],
> > PlotJoined -> True, PlotStyle -> #2,
> > PlotRange -> {{-20, 20}, {0, .2}}]] & @@@ {{data,
> > {}}, {data2,
> > Hue[.7]}}] // Show
> >
> >
> > --
> > Hartmut
> >
>
> There are other flaws in this approach, observe:
>
> If we define the random generator as I did above (following
> suggestions of Bob and others)...
>
> data2 = Table[
> Random[NormalDistribution[0,
> Random[NormalDistribution[5, Sqrt[10]]]]], {10001}];
>
> ...it's (only now, to me) clear that we will have events
> corresponding to a distribution with negative variance. Which
> is nonsense (of course!)
>
> Mathematica complains accordingly:
>
> In[74]:= PDF[NormalDistribution[0, -.5], x]
> From In[74]:=
> NormalDistribution::"posscale": "The scale parameter -0.5 is
> expected \
> to be positive."
> Out[74]=
> PDF[NormalDistribution[0, -0.5], x]
>
> ...but not with Random:
>
> In[119]:= And@@(NumericQ/@data2)
> Out[119]= True
>
> A check like...
>
> In[64]:= data0 = Table[Random[NormalDistribution[0, -5]], {10001}];
> In[66]:= data0plus = Table[Random[NormalDistribution[0, 5]], {10001}];
> In[71]:=
> Block[{$DisplayFunction = Identity, step = 33, scale},
> scale = step/Length[#1];
> ListPlot[{Plus[##]/2, -scale/Subtract[##]} & @@@
> Partition[Take[Sort[#1], {1, -1, step}], 2, 1],
> PlotJoined -> True, PlotStyle -> #2,
> PlotRange -> {{-20, 20}, {0, .2}}]] & @@@
> {{data0plus, {}}, {data0,
> Hue[0]}} // Show
>
> ...makes me assume that
>
> Random[NormalDistribution[0, -a]]
>
> is quite the same as
>
> Random[NormalDistribution[0, a]
>
>
> So we got something, but of doubtful value.
>
>
>
> What is wrong is of course the idea of the variance to be
> itself a normal distributed random variable. We have to
> improve on this. The standard one-sided distribution is the
> Gamma distribution, and we have (e.g.)
>
> In[43]:= Variance[GammaDistribution[5/2, 2]]
> Out[43]= 10
>
> In[44]:= Mean[GammaDistribution[5/2, 2]]
> Out[44]= 5
>
> or Plot
>
> In[42]:=
> Plot[PDF[GammaDistribution[5/2, 2], x], {x, 0, 20}, PlotRange -> All]
>
> Now we can do what we wanted to achieve in first place
>
> In[45]:=
> data3 = Table[
> Random[NormalDistribution[0,
> Random[GammaDistribution[5/2, 2]]]], {10001}];
>
> In[107]:=
> g3 = Block[{step = 50, scale}, scale = step/Length[#1];
> ListPlot[{Plus[##]/2, -scale/Subtract[##]} & @@@
> Partition[Take[Sort[#1], {1, -1, step}], 2, 1],
> PlotJoined -> True, PlotStyle -> #2,
> PlotRange -> {{-20, 20}, {0, .2}}]] & @@@ {{data3, Hue[.7]}}
>
> The distribution of data3 is quite similar to data2, but
> somewhat less accentuated near 0 (of course!).
>
> However we now are in a position to use Mathematica more
> intelligently:
> we may integrate for the distribution searched for
>
> In[97]:= dist[x_]=Integrate[
> PDF[NormalDistribution[0,y],x]*
> PDF[GammaDistribution[5/2,2],y],{y,-\[Infinity],\[Infinity]},
> Assumptions\[Rule]{x\[Element]Reals,x^2>0}]
> Out[97]=
> (1/(6*Pi))*(Sqrt[2*Pi]*HypergeometricPFQ[{},
> {-(1/4), 1/4}, -(x^2/32)] +
> (1/(2^(3/4)*(1/x^2)^(3/4)))*((1/2 + I/2)*
> Gamma[-(3/4)]*HypergeometricPFQ[{}, {1/2, 7/4},
> -(x^2/32)]) - (1/(2^(1/4)*(1/x^2)^(5/4)))*
> ((1/8 - I/8)*Gamma[-(5/4)]*HypergeometricPFQ[{},
> {3/2, 9/4}, -(x^2/32)]))
>
> The expression is complex, certainly because of issues of
> branch cuts and effective integration paths. I cannot delve
> into a discussion of these points, but observe that taking
> the real part
>
> In[111]:=
> g = Plot[N[Re[dist[x]]], {x, -20, 20}, PlotStyle -> {Hue[0]}]
>
> In[112]:= Show[g3, g]
>
> fits well to the distribution as randomly generated.
>
> --
> hw
>
Another note to add: integrating over {y,-\[Infinity],\[Infinity]} is a
little too audacious and no wonder that problems arise.
If we do
In[6]:=
dist[x_, Plus] =
Integrate[
PDF[NormalDistribution[0, y], x]*PDF[GammaDistribution[5/2, 2], y],
{y, 0, \[Infinity]}, Assumptions -> {x > 0}]
Out[6]=
(1/(6*Pi))*(Sqrt[2*Pi]*HypergeometricPFQ[{},
{-(1/4), 1/4}, -(x^2/32)] - (1/(8*2^(3/4)))*
(x^(3/2)*(-4*Gamma[-(3/4)]*HypergeometricPFQ[{},
{1/2, 7/4}, -(x^2/32)] + Sqrt[2]*x*
Gamma[-(5/4)]*HypergeometricPFQ[{}, {3/2, 9/4},
-(x^2/32)])))
In[7]:=
dist[x_, Minus] =
Integrate[
PDF[NormalDistribution[0, y], x]*PDF[GammaDistribution[5/2, 2], y],
{y, 0, \[Infinity]}, Assumptions -> {x < 0}]
Out[7]=
(1/(6*Pi))*(Sqrt[2*Pi]*HypergeometricPFQ[{},
{-(1/4), 1/4}, -(x^2/32)] -
(1/(15*(-(1/x^2))^(3/4)))*(4*2^(1/4)*Pi*
((5*HypergeometricPFQ[{}, {1/2, 7/4}, -(x^2/32)])/
Gamma[-(1/4)] + (Sqrt[2]*HypergeometricPFQ[{},
{3/2, 9/4}, -(x^2/32)])/(Sqrt[-(1/x^2)]*
Gamma[-(3/4)]))))
The second function dist[x, Minus] for x < 0 is unusable. What appears to be
the right solution is
dist[Abs[x], Plus]
This is not the same as Re[dist[x]] I had conjectured above, but the
numerical differences are minute for x less than 200., say. For large x both
distributions become numerical unstable and diverge. End of my wisdom.
--
Hartmut
Prev by Date:
**RE: Re: help in generating a gaussian random variable**
Next by Date:
**How do I put Arrow Heads on the ends of the Axis of my plots?**
Previous by thread:
**RE: Re: help in generating a gaussian random variable**
Next by thread:
**RE: RE: Re: help in generating a gaussian random variable**
| |