       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]]]], {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]]]], {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:= PDF[NormalDistribution[0, -.5], x]
> From In:=
> NormalDistribution::"posscale": "The scale parameter -0.5 is
> expected \
> to be positive."
> Out=
> PDF[NormalDistribution[0, -0.5], x]
>
> ...but not with Random:
>
> In:= And@@(NumericQ/@data2)
> Out= True
>
> A check like...
>
> In:= data0 = Table[Random[NormalDistribution[0, -5]], {10001}];
> In:= data0plus = Table[Random[NormalDistribution[0, 5]], {10001}];
> In:=
> 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}} // 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.)
>
> Out= 10
>
> Out= 5
>
> or Plot
>
> In:=
> 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:=
> data3 = Table[
>       Random[NormalDistribution[0,
>
> In:=
> 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:= dist[x_]=Integrate[
>     PDF[NormalDistribution[0,y],x]*
>     Assumptions\[Rule]{x\[Element]Reals,x^2>0}]
> Out=
> (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:=
> g = Plot[N[Re[dist[x]]], {x, -20, 20}, PlotStyle -> {Hue}]
>
> In:= 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:=
dist[x_, Plus] =
Integrate[
{y, 0, \[Infinity]}, Assumptions -> {x > 0}]
Out=
(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*x*
Gamma[-(5/4)]*HypergeometricPFQ[{}, {3/2, 9/4},
-(x^2/32)])))

In:=
dist[x_, Minus] =
Integrate[
{y, 0, \[Infinity]}, Assumptions -> {x < 0}]
Out=
(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*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