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