趣味问题：100个人围成一圈相互射击，最后存活一人的概率

zgg___

2. da[s_] := Module[{n = Length[s]}, Append[Rest[s] Range[n - 1], 0] + s Range[n - 1, 0, -1]];
3. rr[s_] := Prepend[Most[s], 0];
4. gg[k_, n_, s0_] := Module[{s = Table[0, {n}], ds, i}, If[k == 0, Return[Append[Table[0, {n}], 1]]];
5. Do[ds = da[s0[[k - i + 1]]] Binomial[n - 1 - k + i, i - 1]; If[i == 1, ds = rr[ds];]; s += ds;, {i, k}];
6. If[k < n, Do[ds = s0[[k - i + 1]] Binomial[n - 1 - k + i, i]; If[i == 0, ds = rr[ds];]; s += ds;, {i, 0, k}]];
7. Append[s, 0]];
8. m = 30;
9. s1 = {{0, 0, 1}, {0, 2, 0}, {1, 0, 0}}; ss = {{1, 0}}; n = 2;
10. Do[n++; s0 = s1; s1 = Table[gg[i, n, s0], {i, 0, n}]; AppendTo[ss, Most[Last[s1]]];, {m - 2}]; // Timing
11. sk = Map[#/Apply[Plus, #] &, ss]; // Timing
12. ans = {0, 1}; Do[AppendTo[ans, ans.sk[[k]]];, {k, m - 1}]; // Timing
13. ans // ListPlot

复制代码

把代码贴出来吧。 的确是分2步，其中sk就是第2步的输入（例如：用{1, 0}表示2人射击后剩余0人和1人的几率, {1/4, 3/4, 0}表示3人射击后剩余0、1、2人的几率），ans的那行去做10层说的工作。ss是用状态数表示的sk（例如：{1, 0}对应2人时共有1种射击方式的分布，{2, 6, 0}对应3人时所有的8种射击方式的分布），看，初始时n=2，只有一个ss为{1, 0}。 对于要计算ss的第1步，主要由gg来实现，gg[k,n]表示有n个人时，其中k个人开枪（其余人不开枪）的状态结果（s0参数只是为了递推存数使的），我们想要的ss就是gg[n,n]。看，初始时s1 = {{0, 0, 1}, {0, 2, 0}, {1, 0, 0}}，分别表示gg[0,2]，gg[1,2]和gg[2,2]。可以看到，4个gg[0-3,3]都是可以用这3个已知的gg[0-2,2]来算出的，然后就递推吧。这个算法求的是精确值，如果改为数值计算或许效率更高些。 计算结果在0.5上下震荡，还是比较有趣的。

KeyTo9_Fans

“计算结果在0.5上下震荡，还是比较有趣的。” 是的，很有趣的结果。 猜想：计算结果不收敛，上极限大于$0.5$，下极限小于$0.5$。

zeroieme

但波动好象在衰减

wayne

我算的结果好像不太一样，不知道是不是我问题理解的有误。 设p(n)表示n个人经过足够多次数的 敲钟互杀(无自杀情况) 之后，幸存1人的概率，则算得 p(3)=3/4 p(4)=32/43 p(5)=15/22 思路如下：

wayne

n个人，一次敲钟后，有m 个人立刻死亡. 如果2<=m<=n-1 ，所有可能的情况有 C(n,m)*(m-1)^m*m^(n-m) 种 如果m=n，即全部死了，这个情况就相当于错位排列的个数了。 编程序如下：

1. f[1]=1;f[2]=0;
2. ans=Table[f[n]=Sum[Binomial[n,m](m-1)^m m^(n-m)*f[n-m],{m,2,n-1}]/(Sum[Binomial[n,m](m-1)^m m^(n-m),{m,2,n-1}]+Subfactorial[n]),{n,3,200}];

复制代码

图像跟6楼 倒是极其相似的：

mathe

n个人，一次敲钟后，有m 个人立刻死亡. 如果2 wayne 发表于 2012-4-3 10:59

你没有考虑包含了大于m个人死去情况

wayne

18# mathe

如果2<=m<=n-1 ，所有可能的情况有 C(n,m)*(m-1)^m*m^(n-m) 种

嗯 ，这个计数有问题

wayne

终于和zgg的答案完全一致了。 n个人敲一次钟，互相射击，共有(n-1)^n种情况，其中死掉m个人的情况，计数个数设为g(n,m),这个可利用容斥原理， $g(n,m) = C_{n}^{m} \sum _{j=0}^m (j-1)^j (-1)^{j+m} j^{n-j} C_{m}^{j}$ 于是程序为：

1. g[n_, m_] := Binomial[n, m] Sum[(-1)^(m + j) Binomial[m, j] j^(n - j) (j - 1)^j, {j, 0, m}];
2. f[1] = 1; f[2] = 0;ans = Table[f[n] = Sum[g[n, m]*f[n - m], {m, 2, n - 1}]/((n - 1)^n), {n, 3, 50}]

复制代码

f(3)=3/4, f(4)=16/27, f(5)=15/32, f(6)=11704/28125

wayne

6# zgg___ 计算式子算是找到了，可似乎并不能明显的提高速度。 很好奇zgg是怎么算到300多的

mathe

n个人一次射击有$(n-1)^n$种等概率情况。 设其中1～k号死亡,k+1~n号存活的方案总数为$u(n,k),k>=2$,于是一轮射击正好死亡k个人的方案总数为$C_n^ku(n,k)$ 由于我们知道所有k+1~n号存活的方案总数为$k^{n-k}(k-1)^k$ 于是我们得出 $k^{n-k}(k-1)^k=\sum_{h=0}^{k-2} C_k^h u(n,k-h)$ 或者说，我们得出递推公式$u(n,k)=k^{n-k}(k-1)^k-\sum_{h=1}^{k-2}C_k^h u(n,k-h)$ 特别的，$u(n,2)=2^{n-2}$