王丙參,魏艷華

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

M-H算法的建議分布選擇與比較

王丙參,魏艷華

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

文章比較研究了M一H算法與舍選法的聯(lián)系,給出了建議分布的選擇標(biāo)準(zhǔn)與幾種建議選擇,并比較了這幾種建議分布對(duì)馬氏鏈的影響,舉例說明了M一H算法在貝葉斯推斷中可行、穩(wěn)定、有效。

M一H算法;舍選法;隨機(jī)數(shù);建議分布;貝葉斯估計(jì)

0 引言

隨機(jī)模擬是研究復(fù)雜系統(tǒng)的常用手段,其難點(diǎn)在于目標(biāo)分布的計(jì)算機(jī)仿真,即產(chǎn)生目標(biāo)分布隨機(jī)數(shù)[1]。比如,產(chǎn)生多變量非標(biāo)準(zhǔn)形式且各分量之間相依分布的隨機(jī)數(shù)往往非常困難,難于實(shí)現(xiàn),另外,在不完全數(shù)據(jù)場(chǎng)合進(jìn)行統(tǒng)計(jì)推斷,往往需要產(chǎn)生給定截?cái)喾植茧S機(jī)數(shù),但截?cái)喾植紖?shù)的滿條件分布往往沒有顯示表達(dá)式,直接抽樣也很難。MCMC方法就是一種解決上述問題的簡單有效的貝葉斯計(jì)算方法[2,3],而M-H算法是MCMC的核心,已被列為影響工程技術(shù)與科學(xué)發(fā)展的十大算法之首,對(duì)其研究具有重要的實(shí)際意義[4]。在M-H算法具體操作中,建議分布選擇是否恰當(dāng)直接關(guān)系到馬氏鏈的性質(zhì),比如混合性、收斂性等,進(jìn)而影響后面各種統(tǒng)計(jì)計(jì)算[5,6]?；诖?本文比較研究了M-H算法與舍選法的聯(lián)系,給出了建議分布的選擇標(biāo)準(zhǔn)與幾種建議選擇,并比較了這幾種建議分布對(duì)馬氏鏈的影響,舉例說明了M-H算法在貝葉斯推斷中可行、穩(wěn)定、有效。

1 對(duì)照舍選抽樣法研究M-H算法的基本思想與步驟

當(dāng)某目標(biāo)密度函數(shù)π(x)可被計(jì)算,但不易直接抽樣時(shí),利用MCMC方法可獲得π(x)的一個(gè)近似樣本,基本思想是[5]:通過構(gòu)造一個(gè)非周期不可約的馬氏鏈{X(t)，t=1，2，…},使其平穩(wěn)分布恰好為目標(biāo)分布π(x),這樣便可生成近似服從π(x)的樣本,從而進(jìn)行各種統(tǒng)計(jì)推斷。M-H算法是重要性抽樣的實(shí)現(xiàn),也是一種非常通用的構(gòu)造馬氏鏈方法,它的動(dòng)機(jī)是推廣舍選抽樣法。在舍選抽樣中,有目標(biāo)分布π(x)、建議概率密度函數(shù)q(x)和一個(gè)接受準(zhǔn)則h(x，y)(概率函數(shù)),同樣本文也假設(shè)π(x)是全樣本空間Ω上的目標(biāo)分布,需要在Ω上產(chǎn)生樣本馬氏鏈{x1，…，xt，…},使它的穩(wěn)定分布恰好為目標(biāo)分布π(x)。首先,本文令q(x，y)表示由狀態(tài)x轉(zhuǎn)移到狀態(tài)y的概率,也常記為q(y|x),稱為建議分布或提案分布,則任意選擇不可約轉(zhuǎn)移概率q(x，y)及函數(shù)0＜α(.，.)≤1,對(duì)任意組合(x，y)，x≠y,則:

形成一個(gè)轉(zhuǎn)移核。注意,建議分布也常常用符號(hào)g(x，y)表示。

在t=0時(shí)刻,取X(0)=x(0),其中x(0)從初始分布h(x)中抽樣且滿足π(x(0))＞0。但是,初始點(diǎn)的選擇會(huì)影響馬氏鏈的性質(zhì),因此為了更快達(dá)到預(yù)期效果,在具體操作時(shí),x(0)一般可取初始估計(jì)值,比如矩估計(jì)等。接著,給定X(t)=x,首先由q(.|x)產(chǎn)生潛在轉(zhuǎn)移y,然后以概率α(x，y)接受y(X(t+1)=y),以概率1-α(x，y)仍停留在x (X(t+1)=x)。

假設(shè)當(dāng)前狀態(tài)X(t)=x(t),M-H算法具體步驟如下:

(1)從建議分布q(.|x(t))生成候選值x*。

(4)令t=t+1,返回第(1)步。

注意,M-H比率R(x(t)，x*)總是有定義,這是因?yàn)棣?x)＞0且q(x，y)＞0。每完成上述一個(gè)循環(huán),則生成一個(gè)目標(biāo)分布隨機(jī)數(shù)。在具體實(shí)現(xiàn)M-H算法時(shí),可選取多個(gè)初始點(diǎn)來檢驗(yàn)得到的輸出結(jié)果是否一致,這也可看作是與最優(yōu)化算法的結(jié)合。顯然,在M-H算法中,建議分布的主要目的是產(chǎn)生狀態(tài)轉(zhuǎn)移,即為每個(gè)狀態(tài)構(gòu)造一個(gè)鄰域,并選中鄰域中某個(gè)鄰居狀態(tài)。要使得M-H算法產(chǎn)生的隨機(jī)數(shù)序列是馬氏鏈,這就需要要求建議分布q(x，y)在整個(gè)樣本空間Ω上有定義,通常,這是很容易實(shí)現(xiàn)的。舍選抽樣對(duì)所有的x有q(.|x)=q(.),是最簡單的情況,但MCMC方法推廣了舍選抽樣法,使建議分布變成了條件概率密度函數(shù),而M-H算法的神奇之處也正在于此。當(dāng)然,為了保證馬氏鏈狀態(tài)產(chǎn)生遍歷性的轉(zhuǎn)移,必須要求在x的某個(gè)鄰域Nx里成立:

M-H算法接下來要做的是,將選中的狀態(tài)與當(dāng)前狀態(tài)相比,按一定的概率接受其中之一作為隨機(jī)序列的下一狀態(tài)。另一個(gè)與舍選抽樣法不同的是,接受概率不但取決于下一步狀態(tài)x(t+1),而且和當(dāng)前狀態(tài)x(t)有關(guān)。

利用因M-H算法產(chǎn)生的隨機(jī)數(shù)序列,保證了X(t+1)只依賴X(t),而和以前產(chǎn)生的隨機(jī)數(shù)無關(guān),所以M-H算法構(gòu)造的鏈滿足馬氏性,即為馬氏鏈,但它是否是非周期不可約的則取決于建議分布q(x，y)的選擇。所謂的不可約是指馬氏鏈從任意狀態(tài)出發(fā)都能以一定的正概率轉(zhuǎn)移到其他狀態(tài),即它的轉(zhuǎn)移是遍歷的。不具有遍歷性意味著馬氏鏈的整個(gè)狀態(tài)空間被分割成幾個(gè)互不相通的子空間,這樣,它的子空間的狀態(tài)就不能迭代到其他子空間中去,進(jìn)而也不可能存在平穩(wěn)分布。如果是非周期不可約的,則M-H算法構(gòu)造的馬氏鏈具有唯一的極限平穩(wěn)分布,且可以證明M-H算法構(gòu)造的馬氏鏈以π(x)為其平穩(wěn)分布。事實(shí)上,M-H算法產(chǎn)生的馬氏鏈?zhǔn)强赡娴?即:

這表示馬氏鏈穩(wěn)定時(shí),由y→x的量等于由x→y的量,這就導(dǎo)致π(x)為其平穩(wěn)分布。

利用M-H算法產(chǎn)生樣本的主要目的是估計(jì)X～π(x)的某一函數(shù)的期望。隨著t的增大,由M-H馬氏鏈產(chǎn)生隨機(jī)變量的分布越來越近似該鏈的平穩(wěn)分布,故利用這種方法可以估計(jì)很多感興趣的量,如期望Eπ(h(X))、方差varπ(h(X))及尾概率P(X＞a),其中a為很大的數(shù)字,由馬氏鏈的極限性質(zhì)可知,這些量的估計(jì)都是強(qiáng)相合的。由于生成的隨機(jī)數(shù)可以重復(fù),并且這些抽樣點(diǎn)出現(xiàn)的頻率可以修正目標(biāo)密度與提案分布密度的差異,所以本文應(yīng)該保留馬氏鏈中的重復(fù)值。

但在實(shí)際應(yīng)用中,我們不知道什么時(shí)候馬氏鏈才收斂到平穩(wěn)分布。在M-H算法中,只有在極限情況下才有X(t)～π(x),但是,對(duì)于任何迭代都是有限的,不可能得到精確的邊際分布,而初始點(diǎn)的選擇對(duì)馬氏鏈?zhǔn)怯杏绊懙?。馬氏鏈還未進(jìn)入穩(wěn)定階段之前的狀態(tài)稱為暫態(tài)階段,也可形象稱為熱機(jī)階段。于是,為了提高估計(jì)精度,我們通常舍去馬氏鏈的前D個(gè)值,即刪除一些初始生成值,也就是所謂的預(yù)燒期,但是關(guān)于預(yù)燒期的選擇沒有固定的標(biāo)準(zhǔn),要根據(jù)具體情況具體分析,可能很大,比如D=20000,也可能很小,比如D=200。

2 建議分布的選擇與比較

通常,我們要求建議分布q(x，y)對(duì)于給定的x容易產(chǎn)生q(x，y)的隨機(jī)數(shù),這樣才能使得抽樣變得簡單方便,否則,這只是將生成目標(biāo)分布π(x)的困難轉(zhuǎn)移到了生成建議分布的困難,并不能解決問題。另外,一個(gè)具有某種特定性質(zhì)的建議分布可以增強(qiáng)M-H算法的效果,并可直接通過接受概率的大小來反映。關(guān)于接受概率,應(yīng)注意:接受概率要合適,具體問題要具體分布,而不是越大越好,因?yàn)樘髸?huì)導(dǎo)致較慢的收斂性。Gelman等建議,當(dāng)參數(shù)為1維時(shí),接受概率應(yīng)略小于0.5,當(dāng)參數(shù)大于5維時(shí),接受概率最好為0.25左右[5-7],當(dāng)然,這只是個(gè)人建議,不同的學(xué)者可能有不同的看法。

不同的形式建議分布q(x，y)產(chǎn)生了相應(yīng)的M-H算法變形,下面給出幾點(diǎn)具體建議供參考。

(1)M選擇:如果對(duì)任意的x，y,建議分布都有q(x，y)=q(y，x),即建議分布滿足對(duì)稱性,這時(shí),M-H算法就是Metropolis算法(M算法),此時(shí),接受概率M算法是MCMC方法的最早起源,由N.Metropolis等人于1953年在研究原子與分子的隨機(jī)運(yùn)動(dòng)問題時(shí)引入的隨機(jī)模擬方法。

特別有,如果q(x，y)=q(|x-y|),則成為隨機(jī)移動(dòng)M算法。對(duì)稱建議分布很常用,例如,當(dāng)給定x,q(x，y)可取正態(tài)分布,均值為x,協(xié)方差陣為常數(shù)陣,一種比較簡單的選擇就是取正態(tài)分布N(x(t)，σ2),關(guān)于σ的選擇,應(yīng)注意:由大樣本性質(zhì),通常后驗(yàn)分布都具有比較好的正態(tài)性,因此選擇正態(tài)分布作為建議分布是合適的,其均值應(yīng)為上個(gè)狀態(tài)的值,而方差的大小決定了馬氏鏈在參數(shù)支撐集上的混合程度,因此σ的選擇決定了建議分布的好壞。當(dāng)增量方差σ太大,大部分候選值被拒絕,這就導(dǎo)致M-H算法的效率太低,而σ太小,大部分候選值被接受,這時(shí)得到的馬氏鏈幾乎就是隨機(jī)游動(dòng),由于馬氏鏈從一個(gè)狀態(tài)到另一個(gè)狀態(tài)跨度太小,從而無法實(shí)現(xiàn)在整個(gè)支撐集上的快速移動(dòng),進(jìn)而導(dǎo)致M-H算法的效率不高?？梢姦姨蠡蛱《紝?dǎo)致M-H算法的效率較低。通常的做法是,在具體抽樣時(shí)監(jiān)視候選值的接受概率,Robert,Gelman等建議最好控制在[0.15，0.5]。

(2)獨(dú)立抽樣:如果建議分布q(x，y)與當(dāng)前狀態(tài)x無關(guān),即q(x，y)=q(y),則由此建議分布導(dǎo)出的M-H算法稱為獨(dú)立抽樣。此時(shí),由建議分布產(chǎn)生一個(gè)獨(dú)立鏈,抽取的每個(gè)候選值都與當(dāng)前值獨(dú)立,M-H比率為并且如果g(x)＞0,π(x)＞0,則得到的馬氏鏈就是非周期不可約的。此處,M-H比率也可以表示重要比率,令其中π(x)是目標(biāo)分布,q(x)為包絡(luò)分布,則這表明w(x(t))遠(yuǎn)大于w(x*)的值時(shí),馬氏鏈將會(huì)長時(shí)間停留在當(dāng)前值上。進(jìn)一步,我們有接受概率因此建議分布q(x)應(yīng)與目標(biāo)分布π(x)相似,且尾部包含π(x),這樣獨(dú)立抽樣的效果才可能好,這也導(dǎo)致獨(dú)立抽樣在實(shí)際中很少單獨(dú)使用。

(3)單變量M-H算法:如果目標(biāo)分布X是高維的,則直接產(chǎn)生整個(gè)X比較困難,這時(shí)可以逐個(gè)抽取各分量,但是這需要用到完全條件分布。令…，Xn),考慮Xi|X～i,i=1，2，…，n的條件分布,選擇適當(dāng)?shù)霓D(zhuǎn)移核q(xi→yi|x～i),對(duì)于給定的X～i=x～i產(chǎn)生下一個(gè)可能的yi,然后以概率:

(4)隨機(jī)游動(dòng)鏈?zhǔn)峭ㄟ^簡單變化M-H算法得到的另一種馬氏鏈,實(shí)現(xiàn)方法如下:首先生成ε～h(ε),其中h為密度函數(shù),令x*=x(t)+ε,就可得一個(gè)隨機(jī)游動(dòng)鏈。在這種情況下,初始分布g(x*|x(t))=h(x*-x(t))。對(duì)于h,一般應(yīng)選擇為正態(tài)分布N(0，σ2)、以原點(diǎn)為球心的球面上的均勻分布、尺度變換后的t分布。正態(tài)分布的方差可以控制候選值的擴(kuò)散程度,如果目標(biāo)分布是雙峰分布時(shí),建議選擇σ可適當(dāng)大點(diǎn),在一般情況下,可令σ=1,即h～N(0，1)可以證明,如果h(x)在0的領(lǐng)域內(nèi)為正且π(x)的支撐集連通,則生成鏈?zhǔn)欠侵芷诓豢杉s的。

如果建議分布不隨時(shí)間t而改變,則稱M-H算法是齊次的。當(dāng)然,本文也可構(gòu)造非齊次的M-H算法,擊跑算法就是非齊次的M-H算法的典型代表,在這種方法中,從x(t)出發(fā)選擇的建議分布由兩步產(chǎn)生:先選擇移動(dòng)方向,然后在此方向上產(chǎn)生移動(dòng)距離。當(dāng)X的狀態(tài)空間非常受限,且其他方法很難尋找所有空間的區(qū)域時(shí),利用擊跑算法可能會(huì)取得比較好的效果。遺憾的是,非齊次的M-H算法收斂性質(zhì)更難確定,甚至無法確定。

3 混合參數(shù)貝葉斯估計(jì)中M-H算法建議分布的比較

假定y為已知的觀測(cè)數(shù)據(jù),待估參數(shù)θ的先驗(yàn)分布為p(θ),則θ的貝葉斯估計(jì)往往依賴后驗(yàn)分布p(θ|y)=cL (θ|y)p(θ),其中常數(shù)c為歸一化常數(shù),未知,L(θ|y)為似然函數(shù)。由于在后驗(yàn)分布p(θ|y)中c未知,所以它不能直接用于統(tǒng)計(jì)推斷,但是本文可利用中M-H算法獲得p(θ|y)的近似樣本,進(jìn)而進(jìn)行各種統(tǒng)計(jì)推斷。在M-H算法中,常常利用先驗(yàn)分布作為建議分布,即建議分布g(θ)=p(θ),并取MH比率由于先驗(yàn)分布的支撐集覆蓋后驗(yàn)分布的支撐集,因此獨(dú)立鏈的平穩(wěn)分布π(θ)即為后驗(yàn)分布p(θ|y)。當(dāng)然,先驗(yàn)分布的選擇會(huì)顯著影響馬氏鏈的表現(xiàn),如果選擇不合適,則會(huì)導(dǎo)致馬氏鏈?zhǔn)諗亢苈?達(dá)不到預(yù)期效果。

假定觀測(cè)數(shù)據(jù)y1，…，y100來自混合正態(tài)分布:

觀測(cè)數(shù)據(jù)的直方圖見圖1,θ的先驗(yàn)分布為U(0，1),本文可利用M-H算法構(gòu)造平穩(wěn)分布等于θ的后驗(yàn)分布的馬氏鏈。已知觀測(cè)數(shù)據(jù)是由θ=0.7利用合成法生成的,故θ后驗(yàn)密度應(yīng)該集中在0.7附近,下面考察建議分布對(duì)馬氏鏈的影響。

圖1 由混合分布θN(6，0.52)+(1-θ)N(9，0.52)生成的100個(gè)觀測(cè)值直方圖

(1)利用Beta(1，1)(即U(0，1))作為建議分布,模擬結(jié)果如圖2。

圖2 不同初始值θ(1)生成的馬氏鏈軌跡(左圖θ(1)=0,右圖θ(1)=1)

圖2 表明馬氏鏈很快離開初始值,且似乎很容易從θ后驗(yàn)支撐集的各個(gè)部分抽取值,這種表現(xiàn)稱為混合性良好。為了減少初始值的影響,本文省略前100次迭代,得到:左圖均值為0.6467,右圖均值為0.6483?？梢?無論初始值為0,還是1,馬氏鏈都收斂,且樣本均值都與真值0.7非常接近,即模擬效果都不錯(cuò)。

(2)利用Beta(2，10)作為建議分布,模擬結(jié)果如圖3。

圖3 不同初始值θ(1)生成的馬氏鏈軌跡(左圖θ(1)=0,右圖θ(1)=1)

圖3 生成的馬氏鏈慢慢離開初始值,由于擺動(dòng)明顯,顯然,此鏈沒有收斂到其平穩(wěn)分布,只能得到難以令人信服的結(jié)果,但是,其后驗(yàn)分布仍是此鏈的極限分布,故長期運(yùn)行此鏈,在原則上是可以估計(jì)θ的后驗(yàn)分布。本文同樣省略前100次迭代,得到:左圖均值為0.5910,右圖均值為0.5966。可見,無論初始值為0,還是1,樣本均值都與真值0.7的差異比較大。

圖2和圖3說明建議分布對(duì)馬氏鏈的性質(zhì)影響顯著,因此M-H算法的重點(diǎn)在于選擇合適的建議分布。

下面考慮用隨機(jī)游動(dòng)鏈生成后驗(yàn)分布。假定候選值θ*=θ(t)+ε,ε～U(-a，a),顯然,在馬氏鏈的迭代過程中,候選值有可能落在區(qū)間[0，1]外,這就導(dǎo)致迭代過程變得沒有實(shí)際意義。在這種情況下,本以為迭代結(jié)果很差,但是總體而言,迭代結(jié)果還是出乎意料的好。在實(shí)際操作中,馬氏鏈?zhǔn)侨绱说纳衿?哪怕候選值不在區(qū)間[0，1],甚至密度函數(shù)值變?yōu)樨?fù)值,迭代結(jié)果還是如此完美,仍在0.7附近波動(dòng),比如,當(dāng)a=1時(shí),迭代結(jié)果如圖4左圖所示,本文嘗試了很多取值,結(jié)果都是如此完美。但是,當(dāng)a=8時(shí),迭代結(jié)果就比較差了,如圖4右圖所示,這時(shí)候選值會(huì)經(jīng)常不在[0，1]區(qū)間,但馬氏鏈基本不會(huì)接受不合理的候選值,因?yàn)椴缓侠淼暮蜻x值會(huì)導(dǎo)致的極大似然函數(shù)值變小。

圖4 θ*=θ(t)+ε且ε～U(-a，a)生成的馬氏鏈,左圖a=1,右圖a=8,初始值為0

解決候選值有可能落在區(qū)間[0，1]外的一種比較粗糙的解決辦法是,當(dāng)θ?[0，1]時(shí),重新令后驗(yàn)值為0.5。一個(gè)更好的辦法是重新參數(shù)化。令現(xiàn)在,關(guān)于U生成一個(gè)隨機(jī)游動(dòng)鏈,候選值u*=u(t)+ε1, ε1～U(-b，b)。下面給出兩種看待重新參數(shù)化的觀點(diǎn)。

(1)在θ空間運(yùn)行馬氏鏈。這時(shí),建議分布g(.|u(t))要通過變換成為θ空間的建議分布。令|J(θ(t))|表示θ→u變換的Jacobi行列式的絕對(duì)值在θ(t)的估值,則候選值θ*的M-H比率為:

(2)在u空間運(yùn)行馬氏鏈。這時(shí),θ的目標(biāo)密度要變成u的目標(biāo)密度,即則候選值U*=u*的M-H比率為:

顯然,在圖5(左)中,馬氏鏈混合性非常好,但是圖5 (右)中,效果就比較差,收斂速度較慢。為了減少初始值的影響,本文省略前100次迭代,得到:左圖均值為0.6476,標(biāo)準(zhǔn)差為0.0492,右圖均值為0.6502,標(biāo)準(zhǔn)差為0.0435,樣本均值都與真值0.7非常接近,右圖波動(dòng)比左圖小,注意,這并不代表右圖的模擬結(jié)果更好,只是表明馬氏鏈混合性不好。

圖5 u空間運(yùn)行馬氏鏈,左圖b=1,右圖b=10,初始值為0

圖6 u空間運(yùn)行馬氏鏈,b=0.01,左圖初始值為0,右圖初始值為1

顯然,當(dāng)b=0.01時(shí),不管初始值為0,還是初始值為1,馬氏鏈運(yùn)行效果都比較差。左圖的均值為0.5987,標(biāo)準(zhǔn)差0.0694,均值明顯偏離真實(shí)值。右圖的均值為0.7083,標(biāo)準(zhǔn)差0.0371,讀者也許很驚訝,均值竟然更接近0.7,認(rèn)為模擬結(jié)果更好,其實(shí)不然,我們不能僅僅從一個(gè)偶然的好結(jié)果就判定整個(gè)方法好,還要考慮結(jié)果的產(chǎn)生過程,顯然,右圖馬氏鏈的收斂性和混合性都很差,模擬效果很差。

在所有方法中,估計(jì)值都在0.65附近,比真值0.7偏小,這主要是由于100個(gè)特定樣本值的性質(zhì)導(dǎo)致的。

在重新參數(shù)化空間生成的隨機(jī)游動(dòng)鏈與原始空間的生成的隨機(jī)游動(dòng)鏈相比,具有很多不一樣的性質(zhì),重新參數(shù)化可以提高馬氏鏈的表現(xiàn)。

[1]王丙參,魏艷華,孫永輝.利用舍選抽樣法生成隨機(jī)數(shù)[J].重慶師范大學(xué)學(xué)報(bào):自然科學(xué)版,2013,(6).

[2]魏艷華,王丙參,孫永輝.分組數(shù)據(jù)場(chǎng)合逆威布爾分布參數(shù)貝葉斯估計(jì)的混合Gibbs算法[J].四川大學(xué)學(xué)報(bào):(自然科學(xué)版),2014,(1).

[3]魏艷華,王丙參.基于混合Gibbs算法定數(shù)截尾逆威布爾分布參數(shù)貝葉斯估計(jì)[J].統(tǒng)計(jì)與決策,2014,(11).

[4]肖柳青,周石鵬.隨機(jī)模擬方法與應(yīng)用[M].北京:北京大學(xué)出版社, 2014.

[5]Givens GH,Hoeting JA.計(jì)算統(tǒng)計(jì)[M].王兆軍,劉民千,鄒長亮等譯.北京:人民郵電出版社,2009.

[6]Robert C P,Casella G.Monte Carlo Statistical Methods[M].北京:世界圖書出版社,2009.

[7]劉軍.科學(xué)計(jì)算中的蒙特卡羅策略[M].唐年勝,周勇,徐亮譯.北京大學(xué)出版社,2009.

(責(zé)任編輯/亦民)

O212

A

1002-6487(2016)20-0022-04

國家自然科學(xué)基金資助項(xiàng)目(61104045);天水師范學(xué)院中青年教師科研資助項(xiàng)目(TSA1404)

王丙參(1983—),男,河南南陽人,碩士,講師,研究方向:金融數(shù)學(xué)、統(tǒng)計(jì)計(jì)算。