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/