你是对的; 蒙特卡洛是行不通的。（在一个简单的模拟中——也就是说，在没有任何简化的情况下准确地再现问题情况——每次迭代将涉及 900 次跳蚤动作。空单元格比例的粗略估计是 $1/e$，这意味着在 $N$ 这样的迭代之后的 Monte-Carlo 估计约为 $1/N 1/e (1 - 1/e) = 0.2325\ldots /N$。要确定答案到小数点后六位，您需要将其估计为在 5.E-7 内，并且要达到 95+%（例如）的置信度，您必须将该精度大约减半至 2.5E-7。求解 $\sqrt(0.2325/N) \lt 2.5E-7$大约给出 $N > 4E12$。这将是大约 3.6E15 次跳蚤移动，每次需要 CPU 的几个滴答声。有了一个可用的现代 CPU，您将需要一整年的（高效）计算。

就分析解决方案而言，可以进行一些简化。（这些也可以用来缩短蒙特卡洛计算。）空单元的预期数量是所有单元的空概率之和。要找到这一点，您可以计算每个单元格的入住人数的概率分布。这些分布是通过对每个跳蚤的（独立的！）贡献求和获得的。这将您的问题减少为在该网格上的任何给定单元格之间沿 30 x 30 网格查找长度为 50 的路径数（一个是跳蚤的原点，另一个是您要计算概率的单元格跳蚤的占用）。

你能不能遍历每个跳蚤占据细胞的概率。也就是说，跳蚤 k 最初在单元格 (i(k),j(k)) 中的概率为 1。经过 1 次迭代后，他在 4 个相邻单元格中的每个单元格中的概率为 1/4（假设他不在边缘或一个角落）。然后下一次迭代，这些季度中的每一个都被依次“涂抹”。经过 50 次迭代后，您就有了一个跳蚤的占领概率矩阵。重复所有 900 个跳蚤（如果您利用对称性，这将减少近 8 倍）并添加概率（您不需要一次存储所有这些跳蚤，只需存储当前跳蚤的矩阵（嗯，除非您是非常聪明，您可能需要一个额外的工作矩阵）和当前的矩阵总和）。在我看来，有很多方法可以在这里和那里加快速度。

这根本不涉及模拟。但是，它确实涉及大量计算。计算出以高概率给出略好于 6 dp 精度的答案所需的模拟大小并找出哪种方法更快。我预计这种方法会在一定程度上击败模拟。

虽然我不反对whuber指出的精确到小数点后 6 位的蒙特卡罗解决方案实际上不可能（或不切实际），但我认为可以实现精确到六位数的解决方案。

首先，在Glen_b 之后，粒子在静止状态下是可交换的，因此监测不同单元格的占用情况就足够了（如充分），因为这也构成了马尔可夫过程。下一个时间步$t+1$的占用分布由当前时间$t$的占用决定完成。编写转换矩阵 $K$ 绝对是不切实际的，但模拟转换很简单。

其次，正如shabbychef所指出的，可以遵循 450 个奇数（或偶数）方格上的占用过程，当仅考虑偶数次时，它仍保留在奇数方格上，即平方马尔可夫矩阵 $K^2$。

第三，原始问题仅考虑在 $50$ 马尔可夫转换之后的零占用频率 $\hat{p}_0$。鉴于起点对于马尔可夫链 $(\mathbf{X}^{(t)})$ 的平稳概率分布具有非常高的值，并且考虑到对所有单元格的单一平均值的关注，$$\ hat{p}_0=\frac{1}{450}\sum_{i=1}^{450}\mathbb{I}_0(X_i^ {(50)})$$我们可以认为实现时间 $t=50$ 的链 $(\mathbf{X}^{(t)})$ 是来自平稳概率分布的实现。这大大降低了计算成本，因为我们可以直接从这个平稳分布 $\mathbf{\pi}$ 进行模拟，这是一个多项分布，其概率与偶数角上的 2、3 和 4 成正比，其他单元格分别在边缘和内部单元格上。

显然，平稳分布直接提供了预期的空单元数 $$\sum_{i=1}^{450} (1-\pi_i)^{450}$$ 等于 $166.1069$，

pot=rep(c(rep(c(0,1),15),rep(c(1,0),15)),15)*c(2,
    rep(3,28),2,rep(c(3,rep(4,28),3),28),2,rep(3,28),2)
pot=pot/sum(pot)
sum((1-pot)^450)-450
[1] 166.1069

这非常接近 166.11 美元的蒙特卡洛近似值 [基于 10⁸ 模拟，在我的机器上花了 14 小时]。但对于 6 位小数来说还不够接近。

正如whuber所评论的，估计值需要乘以 2 才能正确回答问题，因此最终值为 332.2137，

分析方法可能很乏味，我还没有考虑过错综复杂的地方，但这里有一个您可能想要考虑的方法。由于您对 50 次响铃后空的单元格的预期数量感兴趣，因此您需要在“单元格中的跳蚤数量”而不是跳蚤的位置上定义马尔可夫链（参见 Glen_b 的回答，它模拟了跳蚤作为马尔可夫链。正如安迪在对该答案的评论中指出的那样，该方法可能无法得到您想要的。）

具体来说，让：

$n_{ij}(t)$ 是 $i$ 行和 $j$ 列的单元格中的跳蚤数量。

然后马尔可夫链从以下状态开始：

对于所有 $i$ 和 $j$，$n_{ij}(0) =1$。

由于跳蚤移动到四个相邻单元格之一，因此单元格状态的变化取决于目标单元格中​​有多少跳蚤，四个相邻单元格中有多少跳蚤以及它们将移动到该单元格的概率。使用此观察，您可以将每个单元格的状态转换概率写为该单元格状态和相邻单元格状态的函数。

如果您希望我可以进一步扩展答案，但这以及对马尔可夫链的基本介绍应该可以帮助您入门。