Re: Simulation
- To: mathgroup at smc.vnet.net
- Subject: [mg19400] Re: [mg19353] Simulation
- From: "Tomas Garza" <tgarza at mail.internet.com.mx>
- Date: Mon, 23 Aug 1999 13:57:14 -0400
- Sender: owner-wri-mathgroup at wolfram.com
Roberto Diego [r-diego at develnet.es] wrote:
> I need help on how the expression:
>
> data=Table[NestList[#+0.01*Random[Integer,{-1,1}]&,5,1000]];
>
> could be modified in order to get -1 with probability 1/6, 0 with
> p 4/6 and
> 1 with p 1/6 instead of 1/3,1/3,1/3
>
> Thanks a lot
Hola Roberto!
You need to apply a transformation to the uniform distribution. The idea
goes like this: let edf be the distribution function from which you wish to
obtain samples.In your example, edf is a step function with jumps of heights
1/6, 5/6, 1/6 at the points -1, 0, 1, respectively. The edf maps the sample
space values onto the interval [0, 1], so that you need in effect the
inverse function of the edf. A random number in [0, 1] will produce a
value - through the inverse function - with probability distribution edf.
Then, take
In[1]:=
edf = {{-1, 1/6}, {0, 5/6}, {1, 1}}
and add an extra point to the left with probability zero. Reverse the
coordinates, and produce an interpolating function, g, with
InterpolationOrder -> 0. Applying g to a uniformly distributed random number
in [0, 1] will then give a random number with the prescribed distribution
edf:
In[2]:=
fde = Reverse /@ Prepend[edf, {-2, 0}]
Out[2]=
{{0, -2}, {1/6, -1}, {5/6, 0}, {1, 1})
In[3]:=
g = Interpolation[fde, InterpolationOrder -> 0]
Out[3]=
InterpolatingFunction[{{0, 1}}, "<>"]
You may verify that g is indeed the inverse of edf:
In[4]:=
Plot[g[x], {x, 0, 1}]
Taking a sample of size, say, 100000, and calculating the frequencies of
each value gives (just to check the result):
In[5]:=
muest = Table[g[Random[]], {100000}];
In[6]:=
<< Statistics`DataManipulation`
In[7]:=
Frequencies[muest]
Out[7]=
{{16668, -1}, {66713, 0}, {16619, 1}}
Besides, the method is not too slow:
In[8]:=
data1 = Table[NestList[# + 0.01*g[Random[Integer]] &, 5, 1000]]; // Timing
Out[8]=
{0.05 Second, Null}
Your original random walk took roughly the same amount of time:
In[9]:=
data = Table[
NestList[# + 0.01*Random[Integer, {-1, 1}] &, 5, 1000]]; // Timing
Out[9]=
{0.05 Second, Null}
Tomas Garza
Mexico City