史莱姆自爆问题

知乎链接

草丛里蹦出一只史莱姆，它每天会做以下一件事（等概率）：自爆；什么也不做；分裂成两只；分裂成三只

问题：史莱姆最终全部自爆完的概率是多少？

编程复杂度太大，因为史莱姆的个数是指数上升的。

import numpy as np

N = 10  # 天数
M = 4  # 每天的分裂情况数
C = (M - 1) ** N + 1  # 最多多少只


def go(a, ind, cnt, n, left):
    if ind == M - 1:
        cnt[ind] = n - sum(cnt[:M - 1])
        c = int(np.sum(cnt * np.arange(M)))
        if c < C:
            a[c] += (1 / M) ** n
        return
    for i in range(left + 1):
        cnt[ind] = i
        go(a, ind + 1, cnt, n, left - i)
    cnt[ind] = 0


def get(n):
    # n只变多的概率
    a = np.zeros(C)
    cnt = np.zeros(M)
    go(a, 0, cnt, n, n)
    return a


def main():
    a = np.zeros(C)
    a[1] = 1
    for i in range(N):
        b = np.zeros(C)
        for j in range(C):
            if a[j] > 0:
                b += get(j) * a[j]
        a = b
        print('day', i, a[0])


main()

这道题其实应该用方程思想：设一只史莱姆爆炸的概率是p，$p=\frac{1}{4}+\frac{p}{4}+\frac{p^2}{4}+\frac{p^3}{4}$，p=sqrt(2)-1或者1或者-sqrt(2)-1(舍)

当每次等概率做2件事，保持现状和爆炸时，一定会灭亡。
当每次等概率做3件事，爆炸、保持现状、分裂繁殖时，一定会灭亡。
只有当每次等概率做4件事，爆炸、保持现状、一个变两个、一个变三个时，才有可能避免灭亡。

实际上，还有一种估算的方式
进行模拟实验case_count次，每个case迭代turn_count轮，看看最终extinct的case的个数。结果还算精确。

import random

extinct = 0
case_count = 10000
turn_count = 1000
for i in range(case_count):
    cnt = 1
    died = False
    for turn in range(10):
        nex = cnt
        for j in range(cnt):
            inc = random.randint(0, 3)
            nex += inc - 1
        if nex == 0:
            died = True
            break
        else:
            cnt = nex
    if died:
        extinct += 1

print(extinct / case_count)

转载于:https://www.cnblogs.com/weiyinfu/p/9189035.html