[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
Prev by Date:
**Re: Factoring question**
Next by Date:
**RE: Re: help in generating a gaussian random variable**
Previous by thread:
**RE: Re: help in generating a gaussian random variable**
Next by thread:
**RE: Re: help in generating a gaussian random variable**
| |