伽爾頓板實(shí)驗(yàn)小球分布的研究

晉宏?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)從入口處投入一個小球時，小球最后落入哪個狹槽是偶然的；當(dāng)投入大量小球時，可看到最后落入各狹槽的小球數(shù)目不相同，在中央的槽內(nèi)小球數(shù)目最多，離中央越遠(yuǎn)的槽內(nèi)小球越少；當(dāng)小球數(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ī)上對該實(shí)驗(yàn)進(jìn)行了蒙特卡羅模擬，并把模擬得到的小球分布與導(dǎo)出的小球理論分布進(jìn)行了比較.

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

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

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

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

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

在下面的模擬中，我們采用了如下設(shè)置：共有24行釘子（m＝24），奇數(shù)行的釘子數(shù)為21個（n＝11），偶數(shù)行的釘子數(shù)為22個.我們?nèi)×? 000 000個小球（N＝1 000 000），對它們在伽爾頓板中的下落過程進(jìn)行了模擬，得到了小球頻數(shù)按照落點(diǎn)位置分布的統(tǒng)計(jì)直方圖，見圖2.從圖中可看出，在正對小球入口處（中央位置）的小球數(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ī)模擬的小球頻率條形圖的比較

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

把我們導(dǎo)出的小球理論分布的概率密度函數(shù)式（9）兩邊取自然對數(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）.接著使用蒙特卡羅方法對伽爾頓板實(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