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