In Riddler Jenga, blocks are randomly placed on top of the stack until there is a collapse. All the blocks have the same alignment (e.g., east-west), and when a block, its center is picked randomly along the block directly beneath it.
On average, how many blocks must you place so that your tower collapses — that is, until at least one block falls off?
Let's assume that each block has length $1,$ so that if $X_i, \, i \geq 1,$ is the center of the $i$th block, where without loss of generality we assume $X_1=0.$ Since each block's center is uniformly distributed over the length of the block directly beneath it, that is, $X_{i+1} \mid X_i \sim U(X_i-0.5, X_i+0.5).$
After block is placed, we can determine whether or not the tower will fall via a center of mass calculation. If the first $n$ blocks are placed at $X_1=0, X_2, \dots, X_n,$ then the tower will fall if there is some $k = 1, \dots, n-1$ such that $$\left|\frac{X_n + \sum_{i=k+1}^{n-1} X_i}{n-k} - X_k \right| > \frac{1}{2}.$$
We could start off by calculating the probability of stack collapse at particular levels. For instance, \begin{align*}\mathbb{P} \{ \text{collapse at level $3$} \} &= \mathbb{P} \left\{ \left| \frac{X_2 + X_2}{2} \right| > \frac{1}{2} \right\}\\ &= \int_{-0.5}^{0.5} \int_{x_2-0.5}^{x_2+0.5} I( x_2 + x_3 \lt -1 ) + I(x_2 + x_3 \gt 1) \,dx_3 \, dx_2 \\ &= \int_{-0.5}^{0.5} (2x_2 - 0.5)_+ + (0.5 - 2x_2)_+ \,dx_2 = \frac{1}{8} \end{align*} However, these integrals get fairly intense quickly.
EDIT: My initial post and Python script did not properly account for the 0-indexing within Python and so was off by 1. It definitely makes me feel for Michael Bolton (no, not that no talent @$$-clown, ... the one in Office Space). I have corrected the details and Python below.
Rather, we can easily simulate the block stacks up to some arbitrarily large level, then check the condition above at each level, finding the first such level where the condition fails and denoting that as the maximum height. I encoded a small Jupyter notebook with this simulation design, displayed below. Using $S=500,000$ simulations gives an estimated mean of $7.1179$ with $95\%$ percentile confidence interval of $(7.1074,7.1283).$ One welcome tidbit from the histogram is that as expected there are about $62,364$ of the $500,000$ simulations with maximum height of $3$, which is $12.4728\%$ or roughly the analytical $1/8$th probability.
import numpy as np
from matplotlib import pyplot as plt
%matplotlib inline
max_blocks = 50
num_samples = 500000
def stack_falls(X, level):
diffs = [np.abs((X[level] + X[(k+1):level].sum()) /(level - k) - X[k]) for k in range(level)]
return (max(diffs) > 0.5)
def get_max_stack(X):
blocks, samples = X.shape
stack_height = np.zeros((samples,1))
for smpl in range(samples):
level = 1
while not stack_falls(X[:,smpl], level):
level += 1
stack_height[smpl] = level
return stack_height+1 # shift index by 1
U = np.random.uniform(low=-0.5, high=0.5, size=(max_blocks, num_samples))
X = np.concatenate([np.zeros((1,num_samples)), np.cumsum(U, axis=0)])
M = get_max_stack(X)
mean = M.mean()
stdev = M.std()
sem = stdev / np.sqrt(num_samples)
max_ = M.max()
print('The average of the max stack height distribution using {} samples is {}.'.format(num_samples, mean))
print('The std. deviation of the max stack height distribution using {} samples is {}.'.format(num_samples, stdev))
print('The SE of the mean of the max stack height distribution using {} samples is {}.'.format(num_samples, sem))
print('95% confidence limit of the mean is ({},{}).'.format(mean-1.96*sem, mean+1.96*sem))
print('The maximum of the max stack height distribution using {} samples is {}.'.format(num_samples, max_))
plt.hist(M, bins=range(int(max_)))
plt.title('Histogram of max_stack_height using {} samples'.format(num_samples))
plt.show()
hist, bins = np.histogram(M, bins=range(int(max_)))
hist, bins
No comments:
Post a Comment