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