3
$\begingroup$

Let $S[16]$ be a binary array i.e, elements of $S$ are 0/1 with elements $S[i]$ are taken uniformly and independently form $\{0,1\}$. Let $k$ be a random element taken uniformly from $\{0,1\}$. I have run the following process.

  1. Take $a,b,c$ randomly and independently from $[1,16]$.
  2. Change $k=k+S[a]$
  3. $S[b]=k+S[c]$.

Run the whole process (1,2,3) until all elements of $S$ are filled with 0. What is expected number of runs for this event?

All additions are binary i.e, in $Z_2$.

  • 0
    Are you implying that this is independent of the initial content of $S$? Or did you forget to specify that? Or are the $S[i]$ also initially uniformly and independently distributed? Also, I presume you mean that $k$ is uniformly distributed over $\{0,1\}$? And that addition is carried out in $\mathbb Z_2$?2011-07-24
  • 0
    $S[i]$ are initially uniformly and independently distributed over $\{0,1\}$. $K$ is also uniformly distributed over $\{0,1\}$. And additions are in $Z_2$2011-07-24
  • 0
    Best to include that information in the question so people don't have to read through the comments to understand the question.2011-07-24
  • 2
    Note that it doesn't make any difference to your process if the array is shuffled after each iteration of the steps. Thus, it should be sufficient to track the number $n$ of elements of $S$ that are 1. In particular, it should be sufficient to analyze the Markov chain with state $(n,k) \in \{0,\ldots,16\} \times \{0,1\}$.2011-07-24

1 Answers 1

7

Note that it doesn't make any difference to your process if the array is shuffled after each iteration of the steps. Thus, it should be sufficient to track the number $n$ of elements of $S$ that are 1.

In particular, it should be sufficient to analyze the Markov chain with state $(n,k) \in \{0,\ldots,16\} \times \{0,1\}$. The transitions of this chain are given by

$$k_{t+1} = \begin{cases} 1-k_t & \text{with probability }\tfrac{n_t}{16} \\ k_t & \text{otherwise} \end{cases}$$

and

$$n_{t+1} = \begin{cases} n_t+1 & \text{with probability }\tfrac{n_t}{16}\left(1-\tfrac{n_t}{16}\right)\text{ if }k_{t+1}=0 \\ n_t-1 & \text{with probability }\tfrac{n_t}{16}\left(1-\tfrac{n_t}{16}\right)\text{ if }k_{t+1}=0 \\ n_t+1 & \text{with probability }\left(1-\tfrac{n_t}{16}\right)^2\text{ if }k_{t+1}=1 \\ n_t-1 & \text{with probability }\left(\tfrac{n_t}{16}\right)^2\text{ if }k_{t+1}=1 \\ n_t & \text{otherwise.} \end{cases}$$

If $n$ hits $0$, the process ends; all we need is the expected time to hit $(0,0)$ or $(0,1)$ from the given initial state ($n \sim \operatorname{Bin}(16, \frac 1 2)$, $k$ uniform). Let me write out the transition matrix and get back on that...


OK, since the transition matrices are $2N \times 2N$, writing down the matrix for $N=16$ here is rather impractical. Let me show instead a cut-down example with $N = 4$.

In general, it's convenient to write the transition matrix as $M = M_n M_k$, where $M_k$ gives the transition probabilities for step 2 and $M_n$ for step 3 in your process. Encoding the state $(n,k)$ as index $i = 2n+k-1$ and working with column probability vectors, the transition matrices for $N=4$ are

$$M_k = \frac 1 4 \begin{bmatrix} 3 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 3 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 2 & 2 & 0 & 0 & 0 & 0 \\ 0 & 0 & 2 & 2 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 3 & 0 & 0 \\ 0 & 0 & 0 & 0 & 3 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 4 \\ 0 & 0 & 0 & 0 & 0 & 0 & 4 & 0 \end{bmatrix} \text{ and } M_n = \frac 1 {16} \begin{bmatrix} 10 & 0 & 4 & 0 & 0 & 0 & 0 & 0 \\ 0 & 6 & 0 & 4 & 0 & 0 & 0 & 0 \\ 3 & 0 & 8 & 0 & 3 & 0 & 0 & 0 \\ 0 & 9 & 0 & 8 & 0 & 9 & 0 & 0 \\ 0 & 0 & 4 & 0 & 10 & 0 & 0 & 0 \\ 0 & 0 & 0 & 4 & 0 & 6 & 0 & 16 \\ 0 & 0 & 0 & 0 & 3 & 0 & 16 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \end{bmatrix},$$

giving

$$M = \frac 1 {64} \begin{bmatrix} 30 & 10 & 8 & 8 & 0 & 0 & 0 & 0 \\ 6 & 18 & 8 & 8 & 0 & 0 & 0 & 0 \\ 9 & 3 & 16 & 16 & 3 & 9 & 0 & 0 \\ 9 & 27 & 16 & 16 & 27 & 9 & 0 & 0 \\ 0 & 0 & 8 & 8 & 10 & 30 & 0 & 0 \\ 0 & 0 & 8 & 8 & 18 & 6 & 64 & 0 \\ 0 & 0 & 0 & 0 & 3 & 9 & 0 & 64 \\ 0 & 0 & 0 & 0 & 3 & 1 & 0 & 0 \end{bmatrix}.$$

To calculate the expected survival time of the process, we need the fundamental matrix $F = I + M + M^2 + \cdots = (I - M)^{-1}$. The elements $F_{ij}$ give the expected number times the process visits state $i$ if started from state $j$. The calculation of $F$ is straightforward, but yields somewhat ungainly fractions even for $N = 4$; rounded to 2 decimal digits,

$$F \approx \begin{bmatrix} 4.96 & 4.00 & 4.48 & 4.48 & 4.48 & 4.48 & 4.48 & 4.48 \\ 2.40 & 4.00 & 3.20 & 3.20 & 3.20 & 3.20 & 3.20 & 3.20 \\ 3.96 & 4.50 & 6.48 & 5.48 & 5.67 & 5.79 & 5.79 & 5.79 \\ 6.12 & 7.50 & 8.56 & 9.56 & 9.37 & 9.25 & 9.25 & 9.25 \\ 3.36 & 4.00 & 5.01 & 5.01 & 6.65 & 6.04 & 6.04 & 6.04 \\ 3.36 & 4.00 & 5.01 & 5.01 & 5.83 & 6.86 & 6.86 & 6.86 \\ 0.84 & 1.00 & 1.25 & 1.25 & 1.54 & 1.64 & 2.64 & 2.64 \\ 0.21 & 0.25 & 0.31 & 0.31 & 0.40 & 0.39 & 0.39 & 1.39 \end{bmatrix}.$$

Summing over the columns of $F$, we get the expected total survival time for each initial state (except, of course, those with $n=0$, for which the survival time is $0$)

$$E = \mathbf 1^T \; F \approx \begin{bmatrix} 25.21 & 29.25 & 34.31 & 34.31 & 37.15 & 37.65 & 38.65 & 39.65 \end{bmatrix},$$

which, multiplied by the initial state distribution vector $P_{2n+k-1} = {N \choose n} \frac 1 {2^{N+1}}$, gives a mean expected survival time of $\frac{7553}{240} \approx 31.47$ steps.


For $N = 16$, the corresponding calculation gives an expected survival time of about $7083$ steps. I can also give the expected survival times for other values of $N$ if you want.

Update: Here are the expected survival times for various $N$ and initial choices of $k$:

$$\begin{array}{r||c|c|c|} N & k_0 = 0 & k_0\text{ random} & k_0 = 1 \\ \hline 1 & 0.500000000 & 0.750000000 & 1.000000000 \\ 2 & 4.375000000 & 4.500000000 & 4.625000000 \\ 3 & 13.74927043 & 14.02225195 & 14.29523346 \\ 4 & 30.87208333 & 31.47083333 & 32.06958333 \\ 5 & 59.12594334 & 60.12456357 & 61.12318380 \\ 6 & 102.9786630 & 104.3672328 & 105.7558025 \\ 7 & 169.3232530 & 171.0531510 & 172.7830490 \\ 8 & 268.6316739 & 270.6468948 & 272.6621157 \\ 9 & 416.5658801 & 418.8161295 & 421.0663789 \\ 10 & 636.3579461 & 638.8020774 & 641.2462088 \\ 11 & 962.3500552 & 964.9556906 & 967.5613261 \\ 12 & 1445.231835 & 1447.973893 & 1450.715950 \\ 13 & 2159.757981 & 2162.617119 & 2165.476257 \\ 14 & 3216.099739 & 3219.061030 & 3222.022322 \\ 15 & 4776.525496 & 4779.577369 & 4782.629242 \\ 16 & 7079.897276 & 7083.030703 & 7086.164130 \\ 17 & 10477.62774 & 10480.83562 & 10484.04350 \\ 18 & 15486.43610 & 15489.71279 & 15492.98948 \\ 19 & 22865.71932 & 22869.06029 & 22872.40127 \\ 20 & 33730.97949 & 33734.38105 & 33737.78262 \\ \end{array}$$

Note that the middle column ($k_0$ random) is simply the average of the $k_0 = 0$ and $k_0 = 1$ columns; I've included it for completeness.

Ps. I ran the calculations up to $N = 40$, and found that for large $N$, $\exp(2.875+0.3804 N)$ gives a pretty good fit for the middle column. (Specifically, it's approximately the logarithmic least squares fit for $N = 35, \ldots, 40$.)


Addendum: For large $N$, there are approximate numerical methods that can be use to avoid inverting $I-M$. Generally, they're based on the fact that you don't really need $F = (I-M)^{-1}$, but just the vector $FP = (I-M)^{-1}P = P + MP + M^2P + \cdots$ (and that, further, you really only care about the sum of the elements of $FP$, not of $FP$ itself).

For example, you can approximate the series $FP = \sum_{m=0}^\infty M^mP$ by starting with the initial state distribution vector $P$ and repeatedly multiplying it with $M$. Of course, since the mean time to absorption is exponential in $N$, and since you have to sum several times that many terms before they become negligible, that naive approach isn't really practical.

However, if you start to compute that sum, you'll eventually notice that, for large enough $m$, $M^{m+1}P \approx \lambda M^mP$, where $\lambda$ is a constant slightly less than $1$. In fact, $\lambda$ will be the dominant eigenvalue of $M$, and the terms $M^mP$ will approach the corresponding eigenvector $\tilde P$. Thus, there will be some $m_0$ such that, for $m > m_0$, the remaining terms can be well approximated as $M^mP \approx \lambda^{m-m_0} M^{m_0}P$, and thus $\sum_{m=m_0}^\infty M^mP \approx (1-\lambda)^{-1} M^{m_0}P$.

Of course, that trick only really helps if the dominant eigenvalue is significantly closer to 1 than the others, but that's quite typical for chains like these, where the absorbing state ($n = 0$) is only reachable through a small number of relatively rare states ($n = 1$). For examples, I calculated the eigenvalues explicitly for $N=32$: the dominant eigenvalue is about 0.9999997076, while the second-largest eigenvalue is only 0.9692222154. Thus, $M^mP$ will converge to the quasi-stationary distribution $\tilde P$ about 10000 times faster than it will decay towards zero.

  • 0
    Thanks. Waiting for the final answer.2011-07-24
  • 0
    Or $\frac{191915616959698822140954370121074028938435014457429042003}{27095127072672200907704180806809457040439579268546560}$ steps, to be precise ;-)2011-07-24
  • 0
    Ilmari: Will you kindly point out me more preciously about the order and entries of $M_n$ and $M_k$? Thanks for your effort.2011-07-24
  • 0
    The entries of $M_n$ and $M_k$ are just the transition probabilities given above. In particular, the diagonal entries of $M_k$ equal $1-\frac{n}{N}$ (where $n = 1+\lfloor\frac{j-1}{2}\rfloor$ and $j$ is the column index), while the off-diagonal entries (one per column, below the diagonal if $j$ is odd, above if $j$ is even) equal $\frac{n}{N}$. The off-diagonal elements of $M_n$, two rows above and below the diagonal, equal $\frac nN(1-\frac nN)$ for odd columns ($k=0$) and $(\frac nN)^2$ and $(1-\frac nN)^2$ respectively for even columns ($k=1$).2011-07-24
  • 0
    ...and the diagonal elements of $M_n$ equal $1-2 \frac nN (1-\frac nN)$ for odd columns and $1-(\frac nN)^2-(1-\frac nN)^2$ for even columns. In particular, note that $M_k$ is a (right) stochastic matrix, i.e. its columns sum to one, while $M_n$ (and thus also $M$) is stochastic except for its first two columns, whose sum is one minus the probability of transitioning to an $n=0$ state.2011-07-24
  • 0
    Just awesome. +100 for such details explanation. What will be the result if we also include $k$ to be even? Thank you very much.2011-07-25
  • 0
    Starting with $k=0$ (and $n \sim \operatorname{Bin}(N,\frac 1 2)$), the expected time to hit $n=0$ is $\frac{74093}{2400} \approx 30.87$ steps for $N = 4$ and $\frac{383661432692753272314314197676082979202272769891885502581}{54190254145344401815408361613618914080879158537093120} \approx 7080$ steps for $N = 16$. Starting with $k = 1$, the corresponding times are $\frac{76967}{2400} \approx 32.07$ for $N=4$ and $\frac{384001035146042016249503282808213136551467287937830665431}{54190254145344401815408361613618914080879158537093120} \approx 7086$ for $N=16$.2011-07-25
  • 0
    Many many thanks. Now I am trying to understand the whole approach with my statistician friend. I am really grateful to you.2011-07-26
  • 0
    There are some good introductions to analyzing Markov chains available on the web. For example, see chapter 11 (and 11.2, which describes exactly the method I used here, in particular) of Grinstead & Snell's _Introduction to Probability_, which is available for free [here](http://www.dartmouth.edu/~chance/teaching_aids/books_articles/probability_book/book.html) (and in particular, chapter 11 alone is [here](http://www.dartmouth.edu/~chance/teaching_aids/books_articles/probability_book/Chapter11.pdf)).2011-07-26
  • 0
    Note, though, a minor notational stumbling block: Grinstead and Snell use row probability vectors and right stochastic matrices, while I used column vectors and left stochastic matrices (I wrote "right" in a previous comment by mistake) in my answer, so you'll need to transpose everything and change the order of multiplication to make the notation match up properly.2011-07-26
  • 0
    Ok. I will start to read the book. Thanks for the information. Two question.2011-07-26
  • 0
    1. For the case $N=1000$, it is hard to compute inverse of matrix of order 2000. So is there any recursion formula or approximation formula for this ? 2. Instead of binary elements, let $S[i]$ and $k$ are filled independently and uniformly from the integer set $[1,1000000]$. Do the process 1,2,3 with addition are in $Z_{32}$. What will be the expected number of runs that all elements of $S[i]$ and $k$ will be divisible by 4, 8, 16,32 etc? I understand that I am bothering you too much and I apologize for this.2011-07-26
  • 0
    @diamond kukur: You should probably make the $\mathbb Z_{32}$ version a separate question. In principle, similar techniques will work for that too, but the larger state space is likely to make things messier.2011-07-26
  • 0
    @Ilmari: Ok, I can understand the problem for $Z_{32}$ version. Thank you very very much for your kind help.2011-07-27
  • 0
    @Ilmari: A very polite and humble request. Will you kindly look at http://math.stackexchange.com/questions/54722/expectation-of-an-event ?2011-07-31