Re: A conditional random number generation problem (please help me!)
- To: mathgroup at smc.vnet.net
- Subject: [mg65672] Re: A conditional random number generation problem (please help me!)
- From: "Ray Koopman" <koopman at sfu.ca>
- Date: Thu, 13 Apr 2006 02:19:49 -0400 (EDT)
- References: <e1ijuv$3oh$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
rjmachado3 wrote:
> I need to know the formula for the random function that return random
> numbers in a range of a and b integers [a,b] but that obey on a custom
> probability (possibly different!) for each integer number on this [a,b]
> range (of course the sum of all integer number probabilities are = 1!).
> Finally, what i want is the general function formula that simulate the
> random behavior (based on a custom probability value for each integer
> number between the [a,b] range. confuse? i hope not! please help me!!!!
>
> what i know so far is that the function formula for generating a "pure"
> random number between [a,b] range is:
>
> rand()*(b-a)+a
>
> where rand() return a random number between 0 and 1.
Here is a generator based on Mark Fisher's MakeDiscreteRandomFunction
[mg53424]. It avoids the built-in Random[], which is known to have
problems.
relprobs is a list of relative probabities, which should be positive
but need not sum to 1.
vals is a list of the values that the random variable may take on.
In your case, use Range[a,b].
If vals is omitted then the integers 1,2,... are used.
In[1]:=
MakeDiscreteRandomFunk[relprobs_List] :=
MakeDiscreteRandomFunk[relprobs, Range@Length@relprobs]
In[2]:=
MakeDiscreteRandomFunk[relprobs_List, vals_List] :=
With[{M2 = 2^27 - 2, m1 = 2^26 - 1., m2 = 2^26 - 2},
With[{fun = Block[{x,p,v},
{p,v} = Transpose@Sort@Transpose@{relprobs,vals};
Function @@ {x, Piecewise @ Transpose @
{Reverse@v, Map[ x >= # &,
Reverse@Most@#*(M2*m1+m2)/Last@#&[FoldList[Plus,0,p]]]} } ] },
Compile[{{n, _Integer}},
Table[fun[Random[Integer,M2]*m1 + Random[Integer,m2]], {n}] ] ]]
In[3]:= getdata = MakeDiscreteRandomFunk[{1,2,3,2,1}];
In[4]:= getdata[20]
Out[4]= {1,5,4,3,2,3,3,2,2,3,5,2,5,3,3,4,2,5,4,4}
Sometimes we don't need the actual observations, just their frequency
distribution. The following generates the frequencies, always taking
vals = Range@Length@probs.
In[5]:=
MakeDiscreteRandomFreqFunk[relprobs_List] :=
With[{M2 = 2^27 - 2, m1 = 2^26 - 1., m2 = 2^26 - 2},
With[{k = Length@relprobs, fun = Block[{x,p,v},
{p,v} = Transpose@Sort@MapIndexed[{#1,#2[[1]]}&,relprobs];
Function @@ {x, Piecewise @ Transpose @
{Reverse@v, Map[ x >= # &,
Reverse@Most@#*(M2*m1+m2)/Last@#&[FoldList[Plus,0,p]]] } } ] },
Compile[{{n, _Integer}}, Module[{i, f = Table[0,{k}]},
Do[i = fun[Random[Integer,M2]*m1 + Random[Integer,m2]];
f[[i]]++, {n}]; f], {{i,_Integer},{f,_Integer,1}} ] ]]
In[6]:= getfreqs = MakeDiscreteRandomFreqFunk[{1,2,3,4}];
In[7]:= getfreqs[1*^5]
Out[7]= {9956,19892,29983,40169}