版權申明：本文為博主窗戶(Colin Cai)原創，歡迎轉帖。如要轉貼，必須註明原文網址

　　http://www.cnblogs.com/Colin-Cai/p/9790468.html 

　　作者：窗戶

　　QQ/微信：6679072

　　E-mail：[email protected] 

　　北師大版九年級上冊第74頁有如下這題：

　　怕圖片看不清楚，我抄一遍如下：

　　將36個球放入標有 1，2，...，12 這 12個號碼的 12 個盒子中，然後擲兩枚質地均勻的骰子，擲得的點數之和是幾，就從幾號盒子中摸出一個球。為了儘快將球模完，你覺得應該怎樣放球？

　　這道題目可謂用意深遠啊，試分析如下。

　　可能的解答？

　　無論如何，我們先得想想題目是什麼意思。所謂質地均勻的骰子，解讀一下，就是每次擲骰子，擲得1-6點中任何一點的概率均為1/6。

　　那麼，同時擲兩枚骰子呢？

　　假設兩枚骰子分別為A、B，那麼一起擲的結果可能如下：

　　A 1點，B 1點

　　A 1點，B 2點

　　...

　　A 6點，B 6點

　　以上一共36種可能，每種可能概率均等，都是1/36

　　於是，我們很容易知道，兩個骰子一起擲得點數之和的概率：

　　2點和12點的概率是1/36

　　3點和11點的概率是2/36(1/18)

　　4

　　5點和9點的概率是4/36(1/9)

　　6點和8點的概率是5/36

　　7點的概率是6/36(1/6)

 　　注：兩個骰子的點數加在一起不可能是1，所以編號為1的盒子是不可能放球的

　　題目實際上考慮的是拿光所有的球，所需要擲骰子的次數的數學期望。而題目是希望找到這個數學期望最少的放法。

　　於是一個“可能”的解答如下：

　　要想更快的拿完，每個盒子的球數應該是 概率 X 球的總數

　　於是，

　　編號為2和12的盒子裡面各放1個球

　　編號為3和11的盒子裡面各放2個球

　　編號為4和10的盒子裡面各放3個球

　　編號為5和9的盒子裡面各放4個球

　　編號為6和8的盒子裡面各放5個球

　　編號為7的盒子裡面放6個球

　　我想，這應該是出題者希望的解答吧，也就是“標準答案”？

　　否則，這個問題，就太複雜了。很可惜，上述放法並不是最好的。

　　蒙特卡洛方法

　　對於一個具體的放法，這個拿完次數的數學期望是多少呢？

　　一開始我在紙上不斷的畫啊，只見一大堆排列組合、無窮級數鋪天蓋地而來。

　　我的個天啊，先再找條路吧。於是我就去選擇蒙特卡洛方法(Monte Carlo Method)。

　　簡單點說，就是用計算機模擬每次擲骰子取球的過程直到取完。實驗反覆做多次，根據大數定理，對於數學期望所在的任意領域，隨著實驗次數的增加，平均擲骰子數量落到這個領域內的概率趨向於1。

　　上面的原理太數學化，自然也超過了初中生的理解範疇。但利用這個原理，我們並不難用任何我們熟悉的語言寫出這個模擬實驗。

　　關鍵就是如何選擇取哪個盒子，本文中我們選擇可以和題目中一樣，使用兩個骰子，每個骰子產生1~6平均分佈，然後加一起。然後這並不具備一般性，對於一般問題我們可以引入輪盤法。

　　我們就以本文為例子，我們如何選擇這12個盒子。

　　假設我們有一個隨機手段，一次可以產生1~36這36個數的其中一個，並且產生每個數的概率都是1/36

　　那麼，我們可以這樣定：

　　如果產生的數是1，則選擇2號盒

　　如果產生的數是2，則選擇12號盒

　　如果產生的數在3~4裡，則選擇3號盒

　　如果產生的數在5~6裡，則選擇11號盒

　　...

　　如果產生的數在31~36裡，則選擇7號盒

　　這樣，正好符合取每個盒子的概率。

　　以上的取法彷彿是一個拉斯維加斯的輪盤，所以叫輪盤法。

　　我們用上面的原理，來做這麼一個簡化的問題：

　　假設有三個盒子，每次選擇1號盒子和2號盒子的概率為1/10，選擇3號盒子的概率是4/5

　　現在，我們來看10個球不同方法選擇完的選取次數的數學期望。

　　程式碼用Python很容易寫出來：

import random
cnt = 0
for i in range(0,10000):
        a = [1,1,8]
        while True:
                cnt += 1
                n = random.randint(1,10)
                k = 0 if n==1 else 1 if n==2 else 2
                if a[k]>0:
                        a[k] -= 1;
                if sum(a)==0:
                        break

print(cnt)

　　上面程式碼就是用蒙特卡洛方法測1號、2號、3號放的球分別為1、1、8，做10000次實驗來統計。

　　按照之前的“解題邏輯”，1、1、8這種放法應該是數學期望最小的。我們就來驗證一下。

　　執行多次，發現每一次輸出的值都在170000左右，那麼我們猜測數學期望應該也在17左右。

　　我們來驗證驗證1號、2號、3號放的球分別為0、0、10的情況，也就是1號、2號盒子都不放球，10個球全部放在概率最高的3號盒子。

　　上述程式碼只需要把第四行後面的陣列改成[0,0,10]即可

　　執行多次，發現每一次輸出的值都在1250000附近，那麼我們猜測數學期望應該也在12.5左右。

　　居然比1、1、8的數學期望要小？

　　再或者真的是小概率事件必然發生？我們看到的是假象？……

　　反覆做過多次實驗，當然應該是真相了。然而蒙特卡洛方法畢竟有概率的成分在裡面，也就是未必絕對靠譜，於是我們還是要深入去解決這個問題。

　　遞迴

　　鑑於蒙特卡洛方法有些“不靠譜”，我們需要精確計算。

　　以簡單情況為例，

　　假設我們現在有三個盒子，1號盒子取到的概率為0.1，2號盒子取到的概率為0.1，3號盒子取到的概率為0.8，

　　現在我們在1號盒子裡放0個球(未放球)，在2號盒子裡放1個球，在3號盒子裡放9個球，

　　我們現在去研究拿完所經歷的“擲骰子”次數的數學期望。

　　我們借用Python的語法，稱這裡的這個數學期望為mean([0.1,0.1,0.8], [0,1,9])

　　這裡，mean函式帶兩個引數，第一個是各個盒子概率的列表，第二個是各個盒子所放球數的列表。

　　我們考慮第一次選擇盒子(擲骰子)，只可能會有以下三種情況：

 　　選擇每個盒子都有個概率，再加上剛剛已經選擇過的這一次，

　　那麼，有

　　mean([0.1,0.1,0.8],[0,1,9]) = mean([0.1,0.1,0.8],[0,1,9]) * 0.1

　　                                             + mean([0.1,0.1,0.8],[0,0,9]) * 0.1

　　                                             + mean([0.1,0.1,0.8],[0,1,8]) * 0.8

　　                                             + 1

　　等式左右都有mean([0.1,0.1,0.8],[0,1,9])，移下項，再除一下，得到：

　　mean([0.1,0.1,0.8],[0,1,9]) = (1 + mean([0.1,0.1,0.8],[0,0,9]) * 0.1 + mean([0.1,0.1,0.8],[0,1,8]) * 0.8) / (1 - 0.1)

　　以上紅色字型部分，是遍歷所有球數不為0的盒子，將這個盒子裡的球減1所得到的問題的數學期望與盒子的概率相乘， 所有這樣的值的累和；

　　以上紅色背景部分，是遍歷所有的球數位0的盒子，將這個盒子取到的概率累和。

　　這樣就得到一個遞迴。

　　另外，也要考慮這個遞迴的邊界。這個倒也容易，也就是當所有盒子都沒球的時候，這個數學期望當然為0。

　　於是就有了以下的程式碼：

def mean(p, n):
        if sum(n)==0:
                return 0
        return (1 + sum(list(map(lambda i:0 if n[i]==0 else \
                mean(p,list(map(lambda j:n[j] if i!=j else n[j]-1,range(0,len(n)))))\
                * p[i],\
                range(0,len(n))))))\
                / (1 - sum(list(map(lambda x,y:x if y==0 else 0,p,n))))

 　　上面似乎像甄嬛體，說人話：

def mean(p, n):
        if sum(n)==0:
                return 0
        f1 = 1
        f2 = 1
        for i in range(len(n)):
                if n[i]!=0:
                        n2 = n.copy()
                        n2[i] -= 1
                        f1 += p[i]*mean(p,n2)
                else:
                        f2 -= p[i]
        return f1/f2

 　　上面是python3，如果是Python2的話，list是沒有copy方法的，需要先匯入copy模組，使用copy.copy來複制list

　　樹遞迴有著太多重複計算，對於36個球，其計算規模何等誇張，顯然是不現實的。

　　為了可以快速的遞迴，就得做一片快取來記錄之前的計算，從而避免重複計算展開，讓計算時間變的可行。考慮到效率，這裡我用C語言寫，只針對本章開頭這個問題，也就是36個球，放在2~12號盒子(1號因為不可能選到，被我排除了)。

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>
typedef union {
        uint64_t i;
        double p;
}data_t;
data_t * data;

double cal_mean(int *n, int sub_pos, int a_cnt, const double *p, const int *seq)
{
        int i, pos = 0;
        double f1 = 1.0, f2 = 1.0;

        if(a_cnt == 0) {
                n[sub_pos]++;
                return 0.0;
        }
        for(i=0;i<11;i++)
                pos += n[i]*seq[i];
        if(data[pos].i) {
                n[sub_pos]++;
                return data[pos].p;
        }
        for(i=0;i<11;i++) {
                if(n[i]) {
                        n[i]--;
                        f1 += cal_mean(n,i,a_cnt-1,p,seq) * p[i];
                } else {
                        f2 -= p[i];
                }
        }
        f1 /= f2;
        data[pos].p = f1;
        n[sub_pos]++;
        return f1;
}

int main(int argc, char**argv)
{
        int i, n[11], seq[11];
        double mean, p[11];

        data = malloc(sizeof(data_t)*atoi(argv[1]));
        if(data == NULL) {
                return 1;
        }
        for(i=0;i<11;i++) {
                p[i] = (6-(i+1)/2)/36.0;
        }

        while(1) {
                if(11 != scanf("%d%d%d%d%d%d%d%d%d%d%d",
                        &n[0],&n[1],&n[2],&n[3],&n[4],&n[5],
                        &n[6],&n[7],&n[8],&n[9],&n[10]))
                        break;
                for(i=0;i<11;i++)
                        printf("%d ",n[i]);
                seq[0] = 1;
                for(i=1;i<11;i++)
                        seq[i] = seq[i-1]*(n[i-1]+1);
                memset(data,0,sizeof(data_t)*seq[10]*(n[10]+1));
                mean = cal_mean(n,0,36,p,seq);
                printf("%.8lf\n", mean);
        }
        free(data);
        return 0;
}

　　這個程式不細講，用到了一點點技巧，比如Python裡寫的時候，遞迴裡用來表示每個盒子球數的列表是複製出來的，但在這個程式裡，表示每個盒子球數的陣列記憶體是共用的。另外一點，為了方便，main函式裡放每個盒子球數的陣列n和每個盒子取到概率的陣列p都是按照從盒子概率從大到小順序的，也就是可以看成順序是7號盒、6號盒、8號盒、5號盒、9號盒、4號盒、10號盒、3號盒、11號盒、2號盒、12號盒。

　　驗證範圍

　　現在，我們有了數學期望的計算方法。就需要對可能的方法進行驗證。

　　根據排列組合知識，利用插板法，36個一樣的球放進11個盒子，所有放法數量應該有

　　這個數量明顯太過於誇張。我們考慮到2號和12號、3號和11號、4號和10號、5號和9號、6號和8號的取到概率是相同的，

　　考慮到對稱性，也就是說，對於這5對盒子，每一對的兩個盒子裡面的球數互換，數學期望都是一樣的。

　　從而我們的驗證範圍可以做一個下降，

 　　只可惜這個數量還是不太現實。

 　　如果對於其中兩個取到概率不相等的盒子A和B，P(A)>P(B)，但A盒球的數量小於B盒球的數量，我們猜測此時取完的數學期望大於A、B兩盒球數互換情況下取完的數學期望。

　　以上命題成立，但證明起來比較複雜，此處略去。　　

　　也就是，對於數學期望最小的情況，球數的大小順序一定和盒子取到概率大小順序一致(相同概率的盒子球數未必相同)。

　　按照上述引理，再根據概率相同的盒子的對稱性，我們可以得到數學期望最小值發生以下的驗證範圍：

　　7號球數 ≥ 6號球數 ≥ 8號球數 ≥ 5號球數 ≥ 9號球數 ≥ 4號球數 ≥ 10號球數 ≥ 3號球數 ≥ 11號球數 ≥ 2號球數 ≥ 12號球數

　　使用遞迴不難用Scheme利用遞迴寫出以下的程式碼列出滿足上述條件的所有7號球數、6號球數、...12號球數：

(define (list-all n pieces max-num)
 (if (= pieces 1)
  (if (> n max-num) '() (list (list n)))
  (apply append
   (map
    (lambda (a) (map (lambda (x) (cons a x)) (list-all (- n a) (- pieces 1) a)))
    (range 0 (+ (min n max-num) 1))))))
(define (pr lst) (if (null? lst) (newline) (begin (display (car lst)) (write-char #\space) (pr (cdr lst)))))
(for-each pr (list-all 36 11 36))

　　Shell

　　計算數學期望的程式和產生驗證範圍的程式都有了，假設編譯後，計算數學期望的程式編譯後叫cal-mean，產生驗證範圍的程式編譯後叫make-range

　　以下shell就可以得到最後的結果

#!/bin/bash
./make-range >input.txt
./cal-mean $(awk '{x=1;for(i=1;i<=NF;i++)x*=$i+1;if(size<x)size=x}END{print size}' input.txt) <input.txt >output.txt
sort -k 12 -n output.txt | head -n 1

　　運行了半個小時，終於得到了最後的結果：

　　8 6 6 4 4 3 3 1 1 0 0 69.56934392

　　前面的8、6、6、4、4、3、3、1、1、0、0就分別是7號盒、6號盒、8號盒、5號盒、9號盒、4號盒、10號盒、3號盒、11號盒、2號盒、12號盒裡的球數，而69.56934392則是取完的擲骰子次數的數學期望，此數學期望是所有情況下最小。從而這種球的方法也就是題目的解答。

　　然而，如此複雜的過程得到的最終結果真的是這道初中數學題的原意？出題的老師真的出對了題目？

 　　完整測試程式，用一個shell程式來表示如下連結檔案：