[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