连续抛掷一枚标准的骰子（即出现6种点数的概率一样，均为1/6）直至点数和超过12。则点数和最有可能是多少？例如，如果点数依次为3、6、1、5，则点数和为15。
解 一种解法是穷举所有点数系列，看看哪个点数和最有可能。但有另一种更巧妙的方法。假如点数系列之和为14，则最后一次抛骰子的点数不可能是1，否则点数和为13（已经超过12）而终止。对于产生14的所有点数系列，如果最后一次的点数均减少1，则可以产生13（因得到14的最后一次不能是1，所以得到14的最后一次骰子点数只能是2、3、4、5、6，那么若最后一次骰子点数改为1，2，3，4，5，这些点数系列之和就会得到13），因此，出现点数和为13的点数系列至少与出现点数和为14的点数系列一样可能。实际上，13的可能性比14还要多，因为得到13的还能以6结尾（如{3，4，6}等）。以此类推，13比15、16等更有可能，因此，点数和最有可能是13。
题目出处[1]中只说明点数和最有可能是13，但并没有给出这种可能性具体是多少，博主也没想出，留待以后解决。通过Matlab仿真，可以得到概率大概是多少：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 |
Total = zeros(1,6); % 记录各点数之和出现的次数，超过12的点数和只能是13~18 iter = 10000; % 生成iter条点数之和超过12的点数数列 for i =1:iter b = 0; % 投骰子的点数之和 while(b<=12) b = b + unidrnd(6); end switch b case 13 Total(1) = Total(1)+1; % 点数之和为13的出现次数加1 case 14 Total(2) = Total(2)+1; % 点数之和为14的出现次数加1 case 15 Total(3) = Total(3)+1; % 点数之和为15的出现次数加1 case 16 Total(4) = Total(4)+1; % 点数之和为16的出现次数加1 case 17 Total(5) = Total(5)+1; % 点数之和为17的出现次数加1 case 18 Total(6) = Total(6)+1; % 点数之和为18的出现次数加1 end end Total_Frequency = Total/sum(Total) |
运行以上代码三次，得到（注意，每次运行结果会略有不同）：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 |
>> dice_roll Total_Frequency = 0.2839 0.2313 0.1938 0.1449 0.0966 0.0495 >> dice_roll Total_Frequency = 0.2791 0.2360 0.1949 0.1457 0.0941 0.0502 >> dice_roll Total_Frequency = 0.2732 0.2308 0.1952 0.1438 0.1075 0.0495 |
可以看到，点数和为13的频率在0.28左右，而点数和为18的频率在0.05左右，因此推断，出现13~18的概率应该分别是某一确定值。
2013.4.16更新：点数和分别是13~18的概率已解出。
思路是，先把问题降一阶：连续抛掷一枚标准的骰子，直至点数和超过6（只能是7，8，9，10，11，12），求出现的概率。然后经过推导，得出求出现7~12的概率解析式，接着进一步泛化，得到求出现6n+1、6n+2、6n+3、6n+4、6n+5、6n+6的通解（n为正整数），如下面的Matlab代码所示：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 |
function X = dice_roll_calc(n) p = 1/6; X = [1/6 1/6 1/6 1/6 1/6 1/6]'; A = [ p p p p p p; 0 p p p p p; 0 0 p p p p; 0 0 0 p p p; 0 0 0 0 p p; 0 0 0 0 0 p]; B = [ 1, 0, 0, 0, 0, 0; p, 1, 0, 0, 0, 0; p+p*p, p, 1, 0, 0, 0; p+2*p^2+p^3, p+p*p, p, 1, 0, 0; p+3*p^2+3*p^3+p^4, p+2*p^2+p^3, p+p*p, p, 1, 0; p+4*p^2+6*p^3+4*p^4+p^5, p+3*p^2+3*p^3+p^4, p+2*p^2+p^3, p+p*p, p, 1]; for i=1:n X = A*B*X; end |
推导过程中，发现了两点规律，其中之一就是上面代码中的B矩阵中元素系数是二项式系数，另一个规律在下面说明。
上面的Matlab代码的功能是，输入参数n，求得出现6n+1、6n+2、6n+3、6n+4、6n+5、6n+6的概率，即连续抛掷一枚标准的骰子直至点数和超过6n。当n=2时，就是原题之解：
|
1 2 3 4 5 6 7 8 9 10 |
>> dice_roll_calc(2) ans = 0.2793 0.2370 0.1923 0.1456 0.0974 0.0485 |
通过对比可以看到，dice_roll.m运行结果与上面接近。
增大n，可以发现另一个规律。
当n=10时：
|
1 2 3 4 5 6 7 8 9 10 |
>> dice_roll_calc(10) ans = 0.2857 0.2381 0.1905 0.1429 0.0952 0.0476 |
当n=100时：
|
1 2 3 4 5 6 7 8 9 10 |
>> dice_roll_calc(100) ans = 0.2857 0.2381 0.1905 0.1429 0.0952 0.0476 |
可以发现，概率值收敛于固定值。
有什么经过多次迭代会收敛于固定值呢？正则链的极限分布π！后来回头仔细看看dice_roll_calc.m的代码，发现的确可以看作是一个马尔可夫过程。dice_roll_calc.m代码中的第19行：X=A*B*X，改写一下，令左X的转置为Y(n+1)，右X的转置为Y(n)，P为(A*B)的转置，则有：
上式可以理解为从状态n到状态n+1的一次转移，Y(n)是状态n时处于各状态的概率，Y(n+1)是状态n+1时处于各状态的概率。因此，经过相当多次转移后，它处于各个状态的概率趋于稳定，且概率稳定值与初始状态无关，稳态概率是以下方程
的唯一解。
[1] http://www.ima.org.uk/_db/_documents/maths_puzzle_diceroll_solution.pdf
[2] 二项式系数：http://baike.baidu.com/view/1923263.htm
Speak Your Mind