Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I am starting with a list of length n that contains all -1's. For instance, when n = 6 the list is {-1, -1, -1, -1, -1, -1}. I pick a random integer from 1 to n and toggle the element by multiplying by -1. In the given list, if I pick the third element, the new list would be {-1, -1, 1, -1, -1, -1}. I continue in this manner, picking an element randomly and toggling it from plus to minus or minus to plus. I wish to know how long it will take for the list to have all 1's. How long will it take to become {1, 1, 1, 1, 1, 1}?

My code is:

toggle[box_, n_] :=
  Block[{r, h = box},
    r = RandomInteger[{1, n}];
    h[[r]] = -h[[r]];
    h]

n = 6;
trials = 40000;
box = ConstantArray[-1, n];
ans = 
  Table[sbox = box; 
        (NestWhileList[toggle[#, n] &, 
                       sbox, 
                       Total[#] != n &] // Length) - 1,
        {trials}];
Mean[ans] // N

The exact answer for this one is 416 / 5. I am getting around 83 for my simulation, but it is deadly slow. Can I speed it up?

share|improve this question
    
I was going to post an answer, but I don't to post an answer with theoretical results when I am not sure what is already known in theory. Anyway this problem is well known and we are essentially after the distribution of a hitting time of a random walk on a hypercube. Because the only quantity that matters is the Total of this list, we could also analyse the discrete time discrete state Ehrenfest process (which is apparently hard to google). I think the distribution of the hitting time may be quite nice, in which a single call to RandomVariate may suffice... –  Jacob Akkerboom Jan 18 at 23:52
1  
I will probably ask a question on math.SE, if I cannot figure out the distribution myself. –  Jacob Akkerboom Jan 18 at 23:53

3 Answers 3

up vote 6 down vote accepted

It is almost always a bad idea to generate random variates one at a time in a time-sensitive construct in MM. This is a rudimentary way of doing a sim. run, with a million variates queued up:

len = 15
box = 2^len - 1
cnt = 1
Catch[Scan[(If[(box = BitXor[box, #]) == 0, Throw[cnt], cnt++]) &, 
   2^RandomInteger[len - 1, 1000000]]] 

Folding it instead of scanning may well be faster.

If you want to solve directly, say for the case of a box of six:

d = DiscreteMarkovProcess[1, {
    {0, 1, 0, 0, 0, 0, 0},
    {1/6, 0, 5/6, 0, 0, 0, 0},
    {0, 2/6, 0, 4/6, 0, 0, 0},
    {0, 0, 3/6, 0, 3/6, 0, 0},
    {0, 0, 0, 4/6, 0, 2/6, 0},
    {0, 0, 0, 0, 5/6, 0, 1/6},
    {0, 0, 0, 0, 0, 0, 1}
   }];

Mean[FirstPassageTimeDistribution[d, 7]]

(* 416/5 *)

For the general solution:

size[n_] := Module[{mp, sa},
   sa = SparseArray[{Band[{2, 1}] -> Range[1/n, (n - 1)/n, 1/n], 
      Band[{2, 3}] -> Range[(n - 1)/n, 1/n, -1/n], {1, 2} -> 
       1, {n + 1, n + 1} -> 1}, {n + 1, n + 1}];
   mp = DiscreteMarkovProcess[1, sa];
   {mp,Mean[FirstPassageTimeDistribution[mp, n + 1]]}];

s=size[20]

(* 3234139734016/2909907 *)

Which takes only a fraction of a second even on the netbook I'm on right now.

With the returned process, you can do all kinds of sorcery, like graph how many bits were 1 on the way to all 1's:

d = s[[1]];
f = RandomFunction[d, {1, 5000}];
ListLinePlot[f]

enter image description here

share|improve this answer

One can use many simultaneous random chains. It considerably speeds up the simulation. Also we can implement own distribution (see Defining Distributional Generators here)

ClearAll[lifetime]
lifetime /: Random`DistributionVector[lifetime[size_, chains_ : 200], n_, prec_] := 
 Module[{m, t, res, pos},
   Flatten[Table[
      m = ConstantArray[-1, size chains];
      t = 0;
      res = ConstantArray[0, chains];
      While[Min[res] == 0, t++;
       pos = Range@chains + chains RandomInteger[size - 1, chains];
       m[[pos]] *= -1;
       res += UnitStep[-res] UnitStep[Total@Partition[m, chains] - size] t;
       ];
      res, {⌈n/chains⌉}]][[1 ;; n]]
   ];

Here size is the size of the list, chains is the number of simultaneous chains (the default value 200 is optimal) and n is the number of generated numbers. Now we can use it as a regular distribution with RandomVariate

N@Mean@RandomVariate[lifetime[6], 100000]

83.4748

Histogram[RandomVariate[lifetime[6], 100000], {0, 150, 2}, AxesOrigin -> 0]

enter image description here

The speed is several times faster then heropup's uncompiled version because I use packed arrays.

One can put the code above to Compile:

lifetimeCFun = Compile[{{size, _Integer}, {chains, _Integer}, {n, _Integer}}, 
   Module[{m, t, res, pos},
    Flatten[Table[
       m = ConstantArray[-1, size chains];
       t = 0;
       res = ConstantArray[0, chains];
       While[Min[res] == 0, t++;
        pos = Range@chains + chains RandomInteger[size - 1, chains];
        m[[pos]] *= -1;
        res += UnitStep[-res] UnitStep[Total@Partition[m, chains] - size] t;
        ];
       res, {⌈n/chains⌉}]][[1 ;; n]]
    ], CompilationTarget -> "C", RuntimeOptions -> "Speed"];

lifetimeC /: Random`DistributionVector[lifetimeC[size_, chains_ : 200], n_, prec_] := 
   lifetimeCFun[size, chains, n];

The speed of my compiled solution is comparable to heropup's compiled solution.

share|improve this answer
    
+1, nicely done. –  rasher Jan 18 at 23:03

On my machine, your algorithm takes about 30 seconds to run 10000 trials. A faster algorithm is

F[n_] := First@NestWhile[{#[[1]] + 1, BitXor[#[[2]], 2^RandomInteger[n - 1]]} &, 
                    {0, 2^n - 1}, #[[2]] != 0 &]

where $n$ is the length of the list. This takes about 12 seconds on my machine. We can make it even faster by using Compile:

G = Compile[{{n, _Integer}}, First@NestWhile[{#[[1]] + 1, BitXor[#[[2]], 
            2^RandomInteger[n - 1]]} &, {0, 2^n - 1}, #[[2]] != 0 &]]

And then the command Table[G[6], {10000}] runs in about 0.68 seconds on my system. But while this is a lot faster than your method, it doesn't scale well with $n$ because the expected number of steps before the process terminates increases very quickly. It is better to compute the expected value of the stopping time using analytic methods, rather than by simulation, although if the goal is to use simulation, this is the first method that came to mind.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.