晉宏?duì)I 劉美云

（榆林學(xué)院能源工程學(xué)院，陜西 榆林 719000）

1 引言

伽爾頓板實(shí)驗(yàn)可以形象地說明大數(shù)目隨機(jī)事件中的統(tǒng)計(jì)規(guī)律，以及統(tǒng)計(jì)規(guī)律中伴隨的漲落現(xiàn)象.伽爾頓板裝置是在一塊豎直木板的上部規(guī)則地釘上許多釘子，木板的下部用豎直隔板隔成許多等寬的狹槽，從板頂漏斗形的入口處可以投入小球，板前覆蓋玻璃，以使小球留在狹槽內(nèi)［1］.實(shí)驗(yàn)表明：當(dāng)從入口處投入一個(gè)小球時(shí)，小球最后落入哪個(gè)狹槽是偶然的；當(dāng)投入大量小球時(shí)，可看到最后落入各狹槽的小球數(shù)目不相同，在中央的槽內(nèi)小球數(shù)目最多，離中央越遠(yuǎn)的槽內(nèi)小球越少；當(dāng)小球數(shù)目較多時(shí)，重復(fù)該實(shí)驗(yàn)，每次得到的小球分布彼此近似地重合［1，2］.伽爾頓板實(shí)驗(yàn)中大量小球的分布服從一定的統(tǒng)計(jì)規(guī)律，近似于正態(tài)分布，但由于伽爾頓板左右側(cè)面的阻擋限制，該分布的范圍與正態(tài)分布有差別，不是從-∞到＋∞.伽爾頓板實(shí)驗(yàn)中小球分布的函數(shù)解析式是什么，教材和其他文獻(xiàn)中沒有給出［1～5］.我們使用最大熵原理研究了伽爾頓板實(shí)驗(yàn)中小球的分布，得到了小球分布的概率密度函數(shù)解析表達(dá)式；我們還在計(jì)算機(jī)上對(duì)該實(shí)驗(yàn)進(jìn)行了蒙特卡羅模擬，并把模擬得到的小球分布與導(dǎo)出的小球理論分布進(jìn)行了比較.

2 伽爾頓板實(shí)驗(yàn)中小球理論分布的推導(dǎo)

最大熵原理是統(tǒng)計(jì)物理中的一個(gè)基本原理，它指出：一個(gè)宏觀系統(tǒng)的信息熵（廣義熵）在一組約束條件下趨于約束極大值.按照此原理，對(duì)于一個(gè)宏觀系統(tǒng)，如果我們選擇合適的約束條件，利用拉格朗日乘子法等方法計(jì)算其信息熵的約束極大值，原則上可以求出該系統(tǒng)的分布［6］.作為自然界的一個(gè)基本規(guī)律，最大熵原理已在很多領(lǐng)域得到廣泛應(yīng)用［6～8］，下面我們使用最大熵原理對(duì)伽爾頓板實(shí)驗(yàn)中小球分布的具體表達(dá)式進(jìn)行推導(dǎo).

圖1 伽爾頓板實(shí)驗(yàn)裝置

伽爾頓板實(shí)驗(yàn)裝置如圖1所示，圖中黑點(diǎn)代表釘子，下面是狹槽.設(shè)入口處相對(duì)于狹槽的高度為h，以板底中心為坐標(biāo)原點(diǎn)，沿板底為x軸建立坐標(biāo)系，見圖1，原點(diǎn)到板底兩端的距離均為L.設(shè)小球落在坐標(biāo)x處的概率密度為f（x），即落在區(qū)間x—x＋dx之間的概率為f（x）dx，由概率歸一化條件，可得

根據(jù)最大熵原理，信息熵S定義為

從入口處投入小球，則小球在下落過程中先后與許多釘子碰撞，最后落入某一狹槽.設(shè)各個(gè)小球落在板底的位置到入口處的距離平方的平均值為C，則有

按照最大熵原理，小球在板底的分布應(yīng)使得信息熵S在約束條件式（1）和式（3）下取得極大值，這類約束極值問題可使用拉格朗日乘子法解決.根據(jù)拉格朗日乘子法，引入函數(shù)：

式中，α是由約束條件式（1）引入的拉格朗日乘子；β是由約束條件式（3）引入的拉格朗日乘子.

由δF［f（x）］＝0，可計(jì)算得信息熵S在約束條件（1）和式（3）下取極大值的概率密度函數(shù)f（x）為

把式（5）代入式（1），得

由式（6）得

與不同x值對(duì)應(yīng)的誤差函數(shù)erf（x）的值可從一般積分表所附的誤差函數(shù)表中直接查出［1］.

把式（7）代回式（5），即可得伽爾頓板實(shí)驗(yàn)中，小球分布的概率密度函數(shù)為

這樣我們便使用最大熵原理推導(dǎo)出了伽爾頓板實(shí)驗(yàn)中小球理論分布的具體函數(shù)表達(dá)式.

3 計(jì)算機(jī)模擬伽爾頓板實(shí)驗(yàn)

計(jì)算機(jī)模擬實(shí)驗(yàn)是科學(xué)研究的重要手段，它可以克服真實(shí)實(shí)驗(yàn)中遇到的許多困難，彌補(bǔ)實(shí)驗(yàn)儀器不足的缺陷［9］.使用計(jì)算機(jī)模擬伽爾頓板實(shí)驗(yàn)可以方便地改變實(shí)驗(yàn)參數(shù)，便于反復(fù)進(jìn)行多次實(shí)驗(yàn)，并快速得到實(shí)驗(yàn)結(jié)果.

我們使用Matlab語言編寫了計(jì)算機(jī)模擬程序，對(duì)伽爾頓板實(shí)驗(yàn)進(jìn)行了蒙特卡羅模擬，模擬的伽爾頓板實(shí)驗(yàn)裝置形狀如圖1所示，釘子的總行數(shù)和每行的釘子個(gè)數(shù)均可調(diào)整.設(shè)共有m行釘子，奇數(shù)行的釘子數(shù)相同為2n-1個(gè)，偶數(shù)行的釘子數(shù)相同為2n個(gè).從上向下統(tǒng)計(jì)行數(shù)，最上邊的釘子為第一行，且第一行中間的那個(gè)釘子正對(duì)入口處，往下每行釘子交錯(cuò)排開.以第一行中間的釘子為坐標(biāo)原點(diǎn)，沿著第一行釘子為x軸，向右為x軸正方向，豎直向下為y軸的正方向，建立坐標(biāo)系.規(guī)定同一行中相鄰兩個(gè)釘子的距離為1，相鄰的兩行距離也是1，奇數(shù)行兩端的釘子與板邊的距離為1，偶數(shù)行兩端的釘子與板邊的距離為0.5；規(guī)定第一行中間釘子的坐標(biāo)為（0，1），從入口處落下的小球第一次都和坐標(biāo)為（0，1）的釘子相碰，即小球落到第一行時(shí)的坐標(biāo)都為（0，1），在隨后的下落過程中每個(gè)小球依次與下面的每行釘子中的一個(gè)釘子相碰，具體和哪個(gè)釘子相碰，由randn函數(shù)生成的一個(gè)正態(tài)分布隨機(jī)數(shù)決定.模擬中小球總數(shù)取為N，定義一個(gè)N行×2列的矩陣，用來存放這N個(gè)小球的位置；矩陣中的每一行代表一個(gè)小球，矩陣的第一列用來存放小球位置的x坐標(biāo)，第二列用來存放小球位置的y坐標(biāo).

現(xiàn)以小球從第一行下落到第二行為例說明一下模擬過程.當(dāng)生成的隨機(jī)數(shù)0≤randn＜n時(shí)，小球下落到第二行時(shí)位于它在第一行位置的右邊，具體下落到哪個(gè)位置是這樣規(guī)定的：當(dāng)0≤randn＜1時(shí)，小球的橫坐標(biāo)加0.5，縱坐標(biāo)加1；當(dāng)1≤randn＜2時(shí)，小球的橫坐標(biāo)加1.5，縱坐標(biāo)加1；……；當(dāng)n-1≤randn＜n時(shí)，小球的橫坐標(biāo)加n-0.5，縱坐標(biāo)加1.當(dāng)生成的隨機(jī)數(shù)-n≤randn＜0時(shí)，小球下落到第二行時(shí)位于它在第一行位置的左邊，具體規(guī)定如下：當(dāng)-1≤randn＜0時(shí)，小球的橫坐標(biāo)加-0.5，縱坐標(biāo)加1；當(dāng)-2≤randn＜-1時(shí)，小球的橫坐標(biāo)加-1.5，縱坐標(biāo)加1；……；當(dāng)-n≤randn＜-n＋1時(shí)，小球的橫坐標(biāo)加-n＋0.5，縱坐標(biāo)加1.當(dāng)生成的隨機(jī)數(shù)n≤randn≤3n時(shí)，小球下落到第二行的位置規(guī)定為：當(dāng)0≤randn-n＜1時(shí)，小球的橫坐標(biāo)為n-0.5，縱坐標(biāo)加1；當(dāng)1≤randn-n＜2時(shí)，小球的橫坐標(biāo)為n-1.5，縱坐標(biāo)加1；……；當(dāng)2n-1≤randn-n≤2n時(shí)，小球的橫坐標(biāo)為-n＋0.5，縱坐標(biāo)加1.當(dāng)生成的隨機(jī)數(shù)-3n≤randn＜-n時(shí)，小球下落到第二行的位置規(guī)定為：當(dāng)-1≤randn＋n＜0時(shí)，小球的橫坐標(biāo)為-n＋0.5，縱坐標(biāo)加1；當(dāng)-2≤randn＋n＜-1時(shí)，小球的橫坐標(biāo)為-n＋1.5，縱坐標(biāo)加1；……；當(dāng)-2n≤randn＋n＜-2n＋1時(shí)，小球的橫坐標(biāo)為n-0.5，縱坐標(biāo)加1.當(dāng)生成的隨機(jī)數(shù)randn＞3n或randn＜-3n時(shí)，拋棄該隨機(jī)數(shù)，令計(jì)算機(jī)再重新生成一個(gè).小球從第二行下落到第三行等繼續(xù)下落的過程依次類推.

在下面的模擬中，我們采用了如下設(shè)置：共有24行釘子（m＝24），奇數(shù)行的釘子數(shù)為21個(gè)（n＝11），偶數(shù)行的釘子數(shù)為22個(gè).我們?nèi)×? 000 000個(gè)小球（N＝1 000 000），對(duì)它們?cè)谫栴D板中的下落過程進(jìn)行了模擬，得到了小球頻數(shù)按照落點(diǎn)位置分布的統(tǒng)計(jì)直方圖，見圖2.從圖中可看出，在正對(duì)小球入口處（中央位置）的小球數(shù)目最多，離中央位置越遠(yuǎn)處的小球數(shù)目越少，分布形狀近似于正態(tài)分布，但與正態(tài)分布有差別之處，正態(tài)分布的范圍是從-∞到＋∞，而伽爾頓板實(shí)驗(yàn)中小球的分布由于受到板左右側(cè)面的阻擋限制，分布范圍是從板的左邊緣到右邊緣.

圖2 小球落點(diǎn)頻數(shù)統(tǒng)計(jì)直方圖

4 理論推導(dǎo)結(jié)果與計(jì)算機(jī)模擬結(jié)果的比較

為了便于比較，我們把上述模擬得到的小球頻數(shù)除以小球總數(shù)，轉(zhuǎn)換為頻率，從而得到了小球頻率按照小球落點(diǎn)位置分布的數(shù)據(jù)，利用這些數(shù)據(jù)作出了小球頻率按照小球落點(diǎn)位置分布的條形圖，見圖3.

圖3 小球概率分布理論曲線與計(jì)算機(jī)模擬的小球頻率條形圖的比較

把頻率取自然對(duì)數(shù)，可得到小球頻率的自然對(duì)數(shù)與小球落點(diǎn)位置間關(guān)系的數(shù)據(jù).然后使用最小二乘法，對(duì)小球頻率的自然對(duì)數(shù)與小球落點(diǎn)位置間關(guān)系的數(shù)據(jù)進(jìn)行二次曲線擬合，擬合得到的二次曲線方程為

把我們導(dǎo)出的小球理論分布的概率密度函數(shù)式（9）兩邊取自然對(duì)數(shù)，得

式（10）與式（11）比較，得到

由上式可計(jì)算得

函數(shù)（14）為我們導(dǎo)出的小球理論分布概率密度函數(shù)的解析表達(dá)式，將此函數(shù)的關(guān)系曲線與小球頻率按照小球落點(diǎn)位置分布的條形圖作在同一張圖上，見圖3.由圖3可見，理論導(dǎo)出的小球分布關(guān)系曲線（細(xì)線）與模擬得到的小球分布頻率的條形圖符合得較好.這說明我們使用最大熵原理導(dǎo)出的伽爾頓板實(shí)驗(yàn)中小球理論分布函數(shù)解析式較好地符合了實(shí)際情況.

5 結(jié)論

本文使用最大熵原理研究了伽爾頓板實(shí)驗(yàn)小球的分布，推導(dǎo)出了小球落點(diǎn)分布的概率密度函數(shù)解析表達(dá)式，見式（9）.接著使用蒙特卡羅方法對(duì)伽爾頓板實(shí)驗(yàn)進(jìn)行了計(jì)算機(jī)模擬.通過把導(dǎo)出的小球理論分布概率密度函數(shù)解析式作成曲線，并把此函數(shù)曲線與小球頻率按落點(diǎn)位置分布的條形圖作在一起進(jìn)行比較，發(fā)現(xiàn)二者符合得較好，這說明導(dǎo)出的小球分布函數(shù)解析式較為成功.

［1］李椿.熱學(xué) 第二版［M］.北京：高等教育出版社，2008

［2］郝志峰，謝國瑞，汪國強(qiáng).概率論與數(shù)理統(tǒng)計(jì) 第二版［M］.北京：高等教育出版社，2009

［3］彭芳麟.伽爾頓板實(shí)驗(yàn)的計(jì)算機(jī)模擬［J］.大學(xué)物理，2005，24（1）：45～49

［4］廖旭，任學(xué)藻.用二項(xiàng)式分布研究伽爾頓板實(shí)驗(yàn)的分布曲線［J］.實(shí)驗(yàn)科學(xué)與技術(shù)，2006，（1）：79～81

［5］聶燕.高爾頓釘板試驗(yàn)的算法實(shí)現(xiàn)及分析［J］.中國民航飛行學(xué)院學(xué)報(bào)，2008，19（3）：62～64

［6］Banavar J R，Maritan A，Volkov I.Applications of the principle of maximum entropy：from physics to ecology［J］.JournalofPhysics：Condensed Matter，2010，22：063101

［7］Plastino A，Curado E M F.Equivalence between maximum entropy principle and enforcing dU＝TdS［J］.PhysicalReviewE，2005，72：047103

［8］Jin H Y，Luo L F，Zhang L R.Using estimative reaction free energy to predict splice sites and their flanking competitors［J］.Gene，2008，424（1-2）：115～120

［9］彭芳麟.計(jì)算物理基礎(chǔ)［M］.北京：高等教育出版社，2010