Polya壇子抽樣模型及其隨機(jī)模擬

王丙參，魏艷華，張藝馨

（天水師范學(xué)院 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，甘肅 天水 741001）

Polya壇子模型是非常重要的概率模型：壇子中有b個(gè)黑球與r個(gè)紅球，每次任取一球，將原球放回后再加入c個(gè)同色球和d個(gè)異色球，它具有重要的應(yīng)用價(jià)值.[1-4]例如，由其衍生出的Polya分布是一種傳染分布或稱(chēng)為概率傳染分布，在氣候統(tǒng)計(jì)中，Polya分布常用來(lái)擬合霧、雷暴等.記ξn為n次中抽到黑球的次數(shù)，董立華在d=0條件下探討了的極限分布，[1]努爾買(mǎi)買(mǎi)提斯吉在c=0,d＞0條件下得出的極限是取值為0.5的退化分布.[3]很多學(xué)者對(duì)Polya壇子模型的理論分析都是限于特定條件下進(jìn)行的，這是因?yàn)閷?duì)一般Polya壇子模型進(jìn)行理論研究非常困難，甚至無(wú)法得出解析結(jié)論.近幾十年來(lái)，鞅論在諸如金融、保險(xiǎn)和可靠性理論等實(shí)際問(wèn)題中得到了廣泛應(yīng)用，是理論探討的一種新方法.[5-6]另外，隨機(jī)模擬方法是研究復(fù)雜系統(tǒng)的一種有效方法，常被比喻為“最后的方法”，因?yàn)樗山鉀Q其他數(shù)值方法難以或不能解決的問(wèn)題.鑒于此，本文系統(tǒng)探討了Polya壇子模型與常見(jiàn)分布的關(guān)系，以便讀者更深入理解與運(yùn)用常見(jiàn)分布，如超幾何分布與Polya分布，利用鞅論與隨機(jī)模擬方法研究了Polya壇子模型，得出一些有趣的結(jié)論，這對(duì)其它概率模型分析提供了思路與方法，最后給出了隨機(jī)模擬的MATLAB代碼供讀者參考.

1 Polya壇子模型與常見(jiàn)分布的關(guān)系

設(shè)壇子中有b個(gè)黑球，r個(gè)紅球，每次隨機(jī)取出一個(gè)球，取出后將原球放回，再加入c個(gè)同色球和d個(gè)異色球.記Bn為第n次取出的是黑球，Rn為第n次取出的是紅球.[7]若連續(xù)從壇子中取出三個(gè)球，其中兩個(gè)紅球、一個(gè)黑球，則有

顯然，以上概率與黑球在第幾次被抽取有關(guān).Polya壇子模型因抽樣結(jié)果與抽樣過(guò)程有關(guān)而難于理論分析.下面在特殊情況下，探討Polya壇子模型的各種變化.

（1）當(dāng)c=-1，d=0時(shí)，前次抽取結(jié)果會(huì)影響后次抽取結(jié)果，但抽取黑球概率不依賴(lài)抽取次序，這也是抽獎(jiǎng)、抓鬮的理論依據(jù)，即不放回抽樣.

（2）當(dāng)c=0，d=0時(shí)，即為放回抽樣，前次抽取結(jié)果不會(huì)影響后次抽取結(jié)果.

（3）當(dāng)c＞0，d=0時(shí)，每次取出球后，會(huì)增加下次取到同色球的概率，即每次發(fā)現(xiàn)一個(gè)傳染病患者，都會(huì)增加以后再傳染的概率，稱(chēng)為傳染病模型.

顯然，利用對(duì)等性有：此結(jié)論也可用數(shù)學(xué)歸納法很容易證明.

（4）當(dāng)c=0，d＞0時(shí)，稱(chēng)為安全模型，也稱(chēng)為Friedman罐子模型.即每當(dāng)事故發(fā)生后，安全工作就該抓緊，下次再發(fā)生事故的概率就會(huì)減少；反之，則否.

當(dāng)d=0時(shí)，假定進(jìn)行n次抽樣，令ξ表示n次中抽到黑球的次數(shù)，則

其中顯然，這構(gòu)成一概率分布，稱(chēng)為Polya分布.當(dāng)c=-1時(shí)，Polya分布ξ就是在產(chǎn)品的質(zhì)量控制中是常用的超幾何分布.[1]

定理1當(dāng)d=0時(shí)，令表示n次中抽到黑球的次數(shù)，則

證明 令表示第i次抽球中抽到黑球的個(gè)

數(shù)，則是同分布于B(1,p)而非獨(dú)立的序列，顯然， Eξi=np， Dξi=npq.于是，

當(dāng) i＜j時(shí)，

令如果對(duì)于保持常數(shù)，則Polya分布收斂于二項(xiàng)分布.如果并令采用數(shù)學(xué)歸納法很容易證明：[1]

其中，表示n次中抽到黑球等于k次的概率.該結(jié)論也可稱(chēng)為Polya分布ξ的負(fù)二項(xiàng)逼近.

現(xiàn)將Polya分布推廣：如果一個(gè)離散型隨機(jī)變量X的分布律為

其中 β,m＞0，則稱(chēng)X服從參數(shù)為 β,m的Polya分布.當(dāng)r≥1時(shí)，令d=βm代表傳染數(shù)量， β代表相對(duì)傳染，則

由 于其中因此推廣Polya分布的分布列可改寫(xiě)為：

進(jìn)一步有，分布函數(shù)

母函數(shù)

矩母函數(shù)

特征函數(shù)

顯然有：

即當(dāng)時(shí)，Polya分布退化為泊松分布，因此Polya分布是泊松分布的連續(xù)修正.[2]

2 Polya壇子模型的鞅分析

定理2在Polya壇子模型中，當(dāng)d=0時(shí)，令Mn表示第n次抽取后壇子中黑球的比例，則{Mn}是一致可積鞅存在，且

證明 令Xn表示第n次抽取后壇子中的黑球數(shù)，則是一個(gè)非齊次的馬爾可夫鏈(MC)，其轉(zhuǎn)移概率為

因?yàn)橹袑?duì) Xn+1有影響的信息都包含在Xn中，所以

即{Mn}是一個(gè)鞅.又因?yàn)槌闪?，所以{Mn}是一致可積鞅.

設(shè)0＜a1＜a2＜1，Mn＜a1且令

表示n次摸球后第一次比例從小于a1到超越a2的時(shí)刻.令Tm=min{T,m}，則對(duì)于m＞n，由停時(shí)定理可知但是

即因?yàn)樯鲜綄?duì)一切m＞n成立，于是有這說(shuō)明至少以概率紅球的比例永遠(yuǎn)不會(huì)超過(guò)a2.同樣的討論可知紅球的比

例從超過(guò)a2到再一次回到a1的最大概率為我們可知從a1出發(fā)超過(guò)a2，再小于a1,…,an有n個(gè)循環(huán)的概率應(yīng)為可見(jiàn)，比例不會(huì)在a1,a2之間無(wú)限次跳躍,由a1,a2的任意性可知存在，記為 M∞.

因?yàn)閧Mn}是一致可積鞅，故

證畢.

設(shè)事件A出現(xiàn)的概率為θ，為估計(jì)θ做了n次獨(dú)立觀察，其中A出現(xiàn)的次數(shù)X服從b(n,θ)，即假如在試驗(yàn)前對(duì)事件A毫無(wú)所知，Bayes建議用均勻分布U(0,1)作為θ的先驗(yàn)分布.由貝葉斯公式可得θ的后驗(yàn)分布

這就是Be(x+1,n-x+1)分布.假如對(duì)壇子中的黑球的比例P一無(wú)所知，只能假定P服從U(0,1)（先驗(yàn)分布），而后驗(yàn)分布M∞是期望為的貝塔分布，進(jìn)一步可得：M∞服從分布.特別當(dāng)b=r=c=1,d=0時(shí)，M∞服從U(0,1)分布.相當(dāng)于對(duì)壇子中的黑球的比例一無(wú)所知且沒(méi)做試驗(yàn).

3 隨機(jī)模擬

對(duì)一般Polya壇子模型進(jìn)行理論分析很困難，但可利用隨機(jī)模擬方法對(duì)其仿真，進(jìn)而從數(shù)值上探討其性質(zhì).隨機(jī)模擬方法提供了對(duì)復(fù)雜隨機(jī)系統(tǒng)進(jìn)行研究的一種思路，但在軟件實(shí)現(xiàn)時(shí)，需要一定的技巧.隨機(jī)模擬的編程實(shí)現(xiàn)是很多研究者的障礙之一，故下面給出Polya壇子模型的隨機(jī)模擬程序供讀者參考.

當(dāng)b=6，,r=4，c=3，d=2，抽樣次數(shù)n=10時(shí)，一共模擬m=100000次，MATLAB程序如下：clear all;c=3;d=2;n=10;m=10^5;%n:每次模擬的抽樣次數(shù)，m:模擬次數(shù)

mp=m1./m %在n次抽樣中黑球出現(xiàn)次數(shù)0到10的概率

m2=hist(bs,11);%利用hist命令統(tǒng)計(jì)黑球出現(xiàn)次數(shù)mp2=m2./m %當(dāng)m較大時(shí)，mp2同mp

bnp=sum(bx)/m %第n次抽到黑球的概率

xb=[0:1:10];Eb=dot(xb,mp);Db=dot(xb.^2,mp)-Eb^2;

[Eb,Db] %由模擬分布列計(jì)算期望與方差

[mean(bs),var(bs)]%由模擬結(jié)果計(jì)算期望與方差

hist(bs) %模擬結(jié)果直方圖

一次運(yùn)行部分結(jié)果如下：

mp=0.0015； 0.0091； 0.0375； 0.0883； 0.1547；0.2065；0.2103；0.1614；0.0909；0.0334；0.0064 bnp=0.5261；ans=5.4840；3.1615

圖1 10次抽樣中黑球出現(xiàn)次數(shù)ξ的10000次模擬結(jié)果的直方圖

這表明，當(dāng)b=6，,r=4，c=3，d=2時(shí)，10次抽樣中黑球出現(xiàn)次數(shù)ξ的分布列為

ξ取值012345678910概率0.00150.00910.03750.08830.15470.20650.21030.16140.09090.03340.0064

且5.4840，3.1615，第10次取出黑球的概率P(B10)=0.5216，100000次 ξ的觀測(cè)值就是向量bs的取值，其直方圖如圖1所示.

上述結(jié)果很難通過(guò)理論分析求得，但通過(guò)隨機(jī)模擬可得出數(shù)值解，非常有實(shí)際意義.

特別有，當(dāng)b=6,r=4，c=3，d=0時(shí)，10次抽樣中黑球出現(xiàn)次數(shù)ξ就是Polya分布，分布列的模擬結(jié)果如表1所示.計(jì)算Polya分布理論值的MATLAB程序如下：

end

lpb %Polya分布列的理論值

表1 Polya分布（b=6,r=4，c=3，n=10）的理論結(jié)果與模擬結(jié)果

第10次取出黑球的概率P(B10)的理論值為0.6，模擬值為0.5991，期望Eξ的理論值為6，模擬值為6.0030，方差 Dξ的理論值為7.3846，模擬值為7.4041.

由以上結(jié)果可知，模擬結(jié)果與真實(shí)值的誤差較小，非常具有參考價(jià)值.隨機(jī)模擬與理論分析的結(jié)果相互佐證，這也從側(cè)面說(shuō)明：模擬程序可行、準(zhǔn)確，值得參考.如果需要考查Polya壇子模型的其它指標(biāo)，只需對(duì)上述程序進(jìn)行簡(jiǎn)單修改就可實(shí)現(xiàn).

[1]董立華.波利亞(polya)分布[J].德州師專(zhuān)學(xué)報(bào),1999,15(2):90-92.

[2]徐曉嶺,王蓉華,顧苑培.Polya分布在氣候統(tǒng)計(jì)中的應(yīng)用[J].數(shù)理統(tǒng)計(jì)與管理,2008,27(2):215-226.

[3]努爾買(mǎi)買(mǎi)提斯吉,楊紀(jì)龍,米輝.關(guān)于罐子模型一個(gè)極限分布的注記[J].南京師范大學(xué)學(xué)報(bào)(工程技術(shù)版),2007,7(2):87-90.

[4]胡學(xué)平,姚勱.一個(gè)Pólya罐子模型的極限定理[J].安慶師范學(xué)院學(xué)報(bào)(自然科學(xué)版),2007,13(2):7-8.

[5]茆詩(shī)松,程依明,濮曉龍,編著.概率論與數(shù)理統(tǒng)計(jì)教程[M].北京:高等教育出版社,2004:40-45.

[6]王丙參,魏艷華,孫永輝.復(fù)合負(fù)二項(xiàng)風(fēng)險(xiǎn)模型的分布函數(shù)[J].統(tǒng)計(jì)與決策,2014,50(2):66-69.

[7]魏艷華,王丙參.蒙特卡洛積分及其改進(jìn)[J].統(tǒng)計(jì)與決策,2017,53(12):71-73.

[8]張波,張景肖.應(yīng)用隨機(jī)過(guò)程[M].北京:清華大學(xué)出版社,2004:130-160.