7
$\begingroup$

Let $A$ be an array of length 1000 with all entries 0. I want to fill up $A$ with ones using the following approach:

At each iteration I take three random integers $(j_1,j_2,j_3)$ from [1,1000] with replacement and do the following:

  1. Set $A[j_1]=1$

    2a. If $A[j_2]=1$ and $A[j_3]=0$, then set $A[j_3]=1$

    2b. if $A[j_2]=0$ and $A[j_3]=1$, then set $A[j_2]=1$

    2c. If $A[j_2]=A[j_3]=0$, do nothing

What is the expected number of such trials to fill up A with all ones?

With replacement means $j_2$ may be equal to $j_3$ or any previous $j_2$ etc.

  • 0
    Or in C: a[j1] = 1; a[j2] = a[j3] = a[j2] | a[j3]; Step 2abc is setting j2 and j3 to the result of an OR operation.2011-07-24

3 Answers 3

3

This does not answer the question but it might be helpful for others who are still intending to answer:

Test[n_] := (iter = 0; cset = ConstantArray[0, n];     While[Count[cset, 1] != n, (iter = iter + 1;       rand = RandomInteger[{1, n}, {3}]; cset[[rand[[1]]]] = 1;       If[cset[[rand[[2]]]] == 1 ||         cset[[rand[[3]]]] == 1, (cset[[rand[[2]]]] = 1;         cset[[rand[[3]]]]  = 1), Null];)]; Return[iter]); 

Using this mathematica code one can generate an instance of the problem (here $n=1000$). I used this as an Monte-Carlo approach to approximate the average number of steps. In the picture below you can see the step counts in the trials ( I simulated 10.000 instanced of the problem ):

enter image description here

If you average the steps you get $n_{average} \approx 2861$ which also fits with trials I did with lower step counts. So we can expect an analytic solution to yield this within a range of a few steps.

  • 1
    @OP: Holy mother of god accept Nate's excellent answer and not my bad one (!)2011-08-01
6

If I am understanding the problem correctly, it appears that we can express everything in terms of the number $X_n$ of 1s in the array at step $n$. This gives a much simpler Markov chain to analyze.

Let $N = 1000$, and suppose at time $n$ we have $k$ 1s in the array. Suppose first that $A[j_1] = 0$ (happens with probability $(N-k)/N$). We set it to 1, gaining one 1 (so now there are $k+1$ 1s). If $A[j_2]$ and $A[j_3]$ are both 1 (prob $((k+1)/N)^2$) or both 0 (prob $((N-(k+1))/N)^2$), no more 1s are gained. If exactly one of $A[j_2], A[j_3]$ is 1 (prob $2(k+1)(N-(k+1))/N^2$), an additional 1 is gained.

Similarly, if $A[j_1] = 1$ (prob $k/N$) then we gain an additional 1 with probability $k(N-k)/N^2$, or no additional 1s with probability $(k^2 + (N-k)^2)/N^2$.

So putting this all together, we have the transition probabilities for $X_n$ are: $\begin{align*} p(k, k) &= \frac{k(k^2 + (N-k)^2)}{N^3} \\ p(k, k+1) &= \frac{k^2(N-k) + 2 (N-k)(k+1)(N-(k+1))}{N^3} \\ p(k, k+2) &= \frac{(N-k)((k+1)^2 + ((N-(k+1))^2)}{N^3} \end{align*} $

Perhaps someone can double-check that these sum to 1, as I haven't the time right now. But this should at least reduce the problem. There might be a slick martingale solution; I'll edit if I think of something.

5

Just to complete Nate Eldredge's answer. As his answer says, it's enough to consider the state at any point as being the number of ones in the array. Let $\mathsf E[k]$ be the expected number of steps needed to fill the array with ones, starting from a state with $k$ ones. Then, calculations essentially as in his answer (in code comments below) give a recurrence relation for $\mathsf E[k]$, which can be solved for instance with the following code (uses Sage; you can put it in a .sage file and run it), giving the answer: $2861.36107443079$

This agrees with user Listing's simulation.

N = 1000 E = [0]*(N+2) def findE(k):     if k>=N: return 0     p1 = k/N          #Pr. that A[j1]=1     p0 = 1-p1         #Pr. that A[j1]=0     p2 = (k+1)/N      #Pr. that A[j2]=0, assuming A[j1] was 0 and was set to 1     p1d = 2*p1*(1-p1) #Pr. that {A[j2],A[j3]} = {1, 0}, assuming A[j1] was 1     p0d = 2*p2*(1-p2) #Pr. that {A[j2],A[j3]} = {1, 0}, assuming A[j1] was 0     #Now,     #E[k] = 1 + p1*(p1d*E[k+1]+(1-p1d)*E[k]) + p0*(p0d*E[k+2]+(1-p0d)*E[k+1])     #so     #E[k]*(1-p1*(1-p1d)) = 1 + p1*p1d*E[k+1] + p0*(p0d*E[k+2]+(1-p0d)*E[k+1])     return (1 + p1*p1d*E[k+1] + p0*(p0d*E[k+2]+(1-p0d)*E[k+1])) / (1-p1*(1-p1d))  for k in range(N+1, -1, -1): E[k] = findE(k) print n(E[0]) 

If the last statement is changed to print E[0] it prints the exact answer as a rational number which in lowest terms is a 3075-digit integer divided by a 3072-digit integer, clearly not worth reporting — and also dashing any hopes that a human of reasonable patience could arrive at the answer by hand.

  • 0
    @user12290: That's the current step. :-) Otherwise all E[k]'s will be 0. The idea is that at any stage, you always execute one step, and then you find the expected number of remaining steps.2011-06-30