[Date Index] [Thread Index] [Author Index]

RE: Re: help in generating a gaussian random variable

• To: mathgroup at smc.vnet.net
• Subject: [mg35535] 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:22 -0400 (EDT)
• Sender: owner-wri-mathgroup at wolfram.com

> -----Original Message-----
> From: Wolf, Hartmut 
> Sent: Wednesday, July 17, 2002 12:46 PM
> Cc: 'Salman Durrani'
> Subject: [mg35535] 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: [mg35535] [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

--------------------------------------------
Hartmut Wolf
T-Systems
Systems Integration
Consultant DA7
Address: Goebelstr. 1-3, Darmstadt, 64293
Postal address: P.O. box 100645, Darmstadt, 64206
Germany
Phone: +49 6151 820 3610
Fax: +49 6151 820 3330
Mobile: +49 175 572 3554
E-Mail: mailto:hartmut.wolf at t-systems.com
Internet: http://www.t-systems.com