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