Re: Re: need little help - no longer!

*To*: mathgroup at smc.vnet.net*Subject*: [mg23589] Re: [mg23548] Re: [mg23530] need little help - no longer!*From*: Andrzej Kozlowski <andrzej at tuins.ac.jp>*Date*: Sun, 21 May 2000 18:12:56 -0400 (EDT)*Sender*: owner-wri-mathgroup at wolfram.com

on 00.5.20 4:10 PM, Wagner Truppel at wtruppel at uci.edu wrote: > "Thank you!" to both Johannes Ludsteck and Kenneth Levasseur, who > supplied me with solutions to my query. > > As Johannes pointed out, the Statistics`MultiDiscreteStatistics` > package has facilities for drawing random numbers from a multinomial > probability distribution. However, I personally prefer Ken's solution: > > RandomIndex[m_] := > (m /. {{a__, b___} :> Length[{a}] /; Plus[a] >= #}) &[Random[]] > > Once again, thank you both. > Wagner > I agree that Ken's solution is very elegant and surprisingly a lot faster than the one derived from the MultiDiscreteDistributions package, but if you mainly care about efficiency than it is still pretty inefficient. The reason is that each time the function RandomIndex is run, first a random number is chosen (using uniform distribution) and then repeated pattern matchings have to be made before an index is chosen (this time with multinomial distribution). This can be considerably speeded up if one first makes a list of indices containing each with required frequency and then just chooses one of them using Random[] (with uniform distribution). Here is one way to do this. As always we have: In[1]:= m = {0.2, 0.4, 0.1, 0.3} Out[1]= {0.2, 0.4, 0.1, 0.3} This makes a list of indices containing each index with required frequency: indices[m_] := Module[{l = Map[Rationalize, m]}, Flatten[Table @@@ Transpose[{Range[Length[m]], Partition[LCM @@ Map[Denominator, l]*l, 1]}]]] In this case the list of indices with frequencies becomes: In[6]:= ind = indices[m] Out[6]= {1, 1, 2, 2, 2, 2, 3, 4, 4, 4} Given such a list we can now define RandomIndex as RandomIndex := ind[[Random[Integer, {1, Length[ind]}]]] Now In[7]:= N[Count[Table[RandomIndex, {1000}], #] & /@ Range[1, 4]/1000] Out[7]= {0.201, 0.406, 0.103, 0.309} In[8]:= Timing[Table[RandomIndex, {1000}]] // First Out[8]= 0.05 Second (On a 233 mghz Mac PowerBook). This is much faster than all the solutions proposed so far except Ted Ersek's one which, however, in the form he sent it, could not be applied to a case with a large number of possible outcomes. (It can of course be modified to work for arbitrary number of outcomes and then, I think, it may well be the fastest one, but I have not tried doing so and testing it). -- Andrzej Kozlowski Toyama International University JAPAN http://platon.c.u-tokyo.ac.jp/andrzej/ http://sigma.tuins.ac.jp/