Expectation of an event

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.

-
Statement very unclear. I think what you're saying is, if $A[r]=0$ then set $A[r]$ to 1; otherwise, if exactly one of $A[s],A[t]$ is 0 then set that one to 1. Is that right? – Gerry Myerson Jun 25 '11 at 3:04
Alternatively, you may be saying this: If $\langle r,s,t \rangle$ is the vector of random integers, do two things: (1) if $A[r] = 0$, then set $A[r]$ to $1$; (2) if exactly one of $A[s]$ and $A[t]$ is $0$, then set that one to $1$. In this interpretation you always do both (1) and (2); in Gerry's you do (2) only if $A[r]$ was already $1$ before the trial. Which of these did you have in mind? – Brian M. Scott Jun 25 '11 at 4:55
Are $j1,j2,j3$ different? – Andrew Jun 25 '11 at 5:35
I take three random integer r,s,t in [1,1000] with replacement. And 1. Set A[r]=1 2. Set A[s]=1 if A[t]=0 3. A[t]=1 if A[s]=0. – user12290 Jun 28 '11 at 5:00
@user12290: Please do not repost questions, as it makes it confusing to have answers and comments in two different places. If you want to change wording, you should edit the existing question instead. I edited this question to replace it with your text from question 48166. – Nate Eldredge Jun 28 '11 at 12:40

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 ):

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.

-
Thanks!${}{}{}{}$ – t.b. Jun 28 '11 at 13:33
No problem, thanks for mentioning the issue. – Listing Jun 28 '11 at 13:47
@OP: Holy mother of god accept Nate's excellent answer and not my bad one (!) – Listing Aug 1 '11 at 10:52

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.

-
I am impressed, what a strong language :-) – Listing Jun 28 '11 at 15:02
@Listing: Yes, indeed. :-) Actually it's more or less just Python, and indeed .sage scripts are translated to Python scripts before being executed. And all this could also be done in "pure" Python (without Sage) using the Fraction module for exact rational arithmetic, but well… had Sage and used it. Of course, surely Mathematica can do it too. – ShreevatsaR Jun 28 '11 at 15:10
Thanks for the solution. It works nicely. Just one quarry. In the recursive formula of $E[k]$, why did you add 1? That is $E[k]=1 + \ldots$. – user12290 Jun 30 '11 at 15:58
@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. – ShreevatsaR Jun 30 '11 at 16:27
Thank you again for the explanation. – user12290 Jul 24 '11 at 11:03

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.

-