AW: Re: Re: need little help - no longer!
- To: mathgroup at smc.vnet.net
- Subject: [mg23605] AW: [mg23589] Re: [mg23548] Re: [mg23530] need little help - no longer!
- From: Wolf Hartmut <hwolf at debis.com>
- Date: Wed, 24 May 2000 02:16:06 -0400 (EDT)
- Sender: owner-wri-mathgroup at wolfram.com
Dear Andrzej,
a remark has to be made, see below:
-----Ursprüngliche Nachricht-----
Von: Andrzej Kozlowski [SMTP:andrzej at tuins.ac.jp]
Gesendet am: Montag, 22. Mai 2000 00:13
An: mathgroup at smc.vnet.net
Betreff: [mg23589] Re: [mg23548] Re: [mg23530] need little
help - no longer!
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/
Certainly I agree with you, the notion of "elegance" is sometimes
misleading. What _is_ elegance after all? Considering a piece of code as
"elegant" if it is pleasant, condensed, while still readable, and containing
an aspect of surprise and wit - so I agree of Ken's code as being elegant.
Yet elegance sometimes may draw one off the problem - as for a chess player
when "playing for the gallery".
Having been off for several day's, I tried on this before peeking to the
answers, and came up with
makeRandomIndex[m_] :=
Module[{m2 = Rest[FoldList[Plus, 0, m]], r, quasiSequence =
Sequence},
With[{code = Which @@ MapThread[
quasiSequence[r <= #1, #2] &, {m2,
Range[Length[m2]]}]},
(r = Random[]; code) &]]
pick = makeRandomIndex[m]
(r$152 = Random[]; Which[r$152 <= 0.2, 1, r$152 <= 0.6, 2, r$152 <=
0.7, 3, r$152 <= 1., 4]) &
tt = Table[pick[], {10000}]; // Timing
{1.592 Second, Null}
N[Count[tt, #] & /@ Range[Length[m]]/10000]
{0.2108, 0.4012, 0.0904, 0.2976}
The point to be respected is the introduction of that ominous
"quasiSequence".
Sequence[r <= #1, #2]& would not do, since that will be interpreted as
Function[r <= #1, #2]. I first had Unevaluated in this place, but this
gives rise to a performance penalty, since Unevaluated remains within the
definition of "pick" and such the internal form of Which has to be (re)built
at runtime for each application. Of course you may introduce Hold or
another head and then replace that with Sequence after Mapping (and before
executing With), yet I consider this to be less "elegant".
On my machine your code was faster (1.102 Second) for this application,
however when I tried it on another example, e.g.
mt = With[{t = Table[Random[], {10}]}, t/Plus @@ t]
{0.0808559, 0.103876, 0.102807, 0.0350509, 0.128651, 0.0261343,
0.139011, 0.0463357, 0.170722, 0.166556}
it refused to work. I isolated two problems:
(1) instead of Rationalize one has to apply Rationalize[#, dx]& but
(2) dx -> Precision[#] again would not work, since Table rejects that
large iterators and even dx = 10^-3 gives rise to gorgeous structures.
So your proposal in fact only works for a smaller class of problems
(rationals with small denominators).
Kind regards, Hartmut