When does a list contain only ones?

Question

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?

Answer 1

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]

Answer 2

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.

Answer 3

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]

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.