第11章-Markov链Monte-Carlo方法.pdf

1/52kJIkJIGoBackFullScreenCloseQuitAAALLL?IIIOOO?2/52kJIkJIGoBackFullScreenCloseQuit11111 MarkovMonteCarlo 11.1 O?Monte Carlo 11.2 MarkovMonte Carlo0 11.3 Metropolis-Hastings 11.4 Gibbs?11.5?“dMCMC?O3/52kJIkJIGoBackFullScreenCloseQuitMonte Carlo=?O3AO!7K!LO!On?+kX2?AIBuffonu1777cJ?|?K|?1?CqL?Vl?wMonte Carlo?d/OI0/I0.John von Neumann?n(5K?Monte Carlo4?u.?x?Monte Carlo54/52kJIkJIGoBackFullScreenCloseQuitMonte Carlo?n?)u)?VC?L?O?,?,Au)?5CqLu)?Vl?K?)3MonteCarlon%KO?EVL!lV?!?O5/52kJIkJIGoBackFullScreenCloseQuitL?E?5O?MonteCarlo?Monte Carlo3p?eO?Monte Carlo?n?Monte Carlo=L?EMarkov?4-5?O?MarkovMonte CarloePMCMCMarkov Chain Monte Carlo?0?O?Monte Carlo9Metropolis-HastingsGibbs?6/52kJIkJIGoBackFullScreenCloseQuit11.1OOO?Monte Carlo?Monte CarloL?E?5OkOf;?Buffon?K 11.1.1 xX?1m?la d/l(l a)?1?V)ZL?:?C1?lL?1?K(,Z)(?k0 Z fraca2,0 =?kU Z?/5Ld/PG7/52kJIkJIGoBackFullScreenCloseQuit?1?ZAvXZ l2sin=11-1(b)?KdPR.11-1 Buffon?1?Vp=L(R)L(G)=12R0lsind12a=2la8/52kJIkJIGoBackFullScreenCloseQuit|d?Y5?NgP?gn2nNVp?u?nN2la=?2lNaneR5yBuffon?1?)Z30,a2S?30,S?=kI)l0,a2!?9l0,!?9/52kJIkJIGoBackFullScreenCloseQuit2?ng?e?1?KL?b?ngK?O2lNanm?la9?l(l touzhen function(N,l=0.75,a=1)+n 0;beta runif(N,0,pi);z runif(N,0,a/2)+for(iin1:N)+if(zi=l sin(betai)/2)+n touzhen(50000)11/52kJIkJIGoBackFullScreenCloseQuit 11.1.2:O?0,1?g(x)3m0,1?R10g(x)dx.)b?CXY l0,1?!pK?(X,Y)l/0,10,1?!Vf(x,y)=(1,0 x 1,0 y?-F/SOBu)?VP(B)=PY g(X)=ZY g(X)f(x,y)dxdy12/52kJIkJIGoBackFullScreenCloseQuitOd?-?P(B)=Z10Zg(x)01dydx=Z10g(x)dx=R10g(x)dx?O8(?Bu)?VdBernoulli-E?Bu)?CqLBu)?V1?)l0,1!?2?ng:?Pw:Y g(X)u)?5CqLu)?V=?R10g(x)dx.XOIO?30,1?12R10ex22dxR5y?L“X13/52kJIkJIGoBackFullScreenCloseQuite mc1 function(n)+k 0;x runif(n);y runif(n)+for(i in 1:n)+if(yi exp(xi2/2)/sqrt(2 pi)+k mc1(10000)14/52kJIkJIGoBackFullScreenCloseQuit,R?SintegrateOd integrate(function(x)exp(x2/2)/sqrt(2 pi),0,1)(11.1.1)0.3413447 with absolute error mc3 function(n)+x runif(n);f 1:n+for(i in 1:n)+fi mc3(10000)O?(J:9?O?C17/52kJIkJIGoBackFullScreenCloseQuit11.2MarkovMonte Carlo0003p?eO?!0?Monte Carlo?n?MonteCarlo=MarkovMonte Carlo?up?MCMCkAMarkov?4-Kl8I?)?l?-G?Markov?)?MarkovATvl?uU?-5MCMC?nA4nkd15n5.3.3H?Markov=18/52kJIkJIGoBackFullScreenCloseQuit?!?!?Markov?4-?-G?mS=1,2,N?MarkovXn,n=0,1,2,=V?PP=p11p12 p1Np21p22 p2N pN1pN2 pNNpij(i,j=1,2,N)L?MarkovlG?iL?=?G?j?VuH?Markov4j=limnpnij,j SjLLm$1?X?uG?j?m?19/52kJIkJIGoBackFullScreenCloseQuit-=(1,2,N)Kj(j=1,2,N)ve?5|?)(1,2,N)=p11p12 p1Np21p22 p2N pN1pN2 pNN(1,2,N)NXj=1j=1P(=PPNj=1j=120/52kJIkJIGoBackFullScreenCloseQuitj(j=1,2,N)?Markov?-V?j=1,2,Nkj=PNi=1ipijLMarkov?uG?j?m?uMarkovlG?i=?G?j?iPNj=1j=1 KLMarkov?uG?j?m?j1.dn?d?$1v?m?Markov?40?)11.2.15?-Vk,?21/52kJIkJIGoBackFullScreenCloseQuitb?3?xj(j=1,2,N)?(xipij=xjpjiPNj=1xj=1kG?i?PNi=1xipij=xjPNi=1pji=xjd11.2.1?)?5 j=xj.?i 6=jkipij=jpjiKdMarkovm_?dLeVj,j=1,2,NLMarkov?Kl?m?X?CzE=Vpij?Markov22/52kJIkJIGoBackFullScreenCloseQuitgMarkovkXe-5nnn 11.2.1?b?Xn,n=0,1,-?H?MarkovzXn?KXn?CX?g?Eg(X)3n kgn=1nnXi=1g(Xi)Eg(X),a.s.n11.2.1Markov?HA?4e?23/52kJIkJIGoBackFullScreenCloseQuitnnn 11.2.2 e?AH=t(0 1)?fu-?Kkq(n)fn Ef(x)N(0,1),n n11.2.2H?Cz?IO?d-?K3$1dv?m?T?-G?d?ul?24/52kJIkJIGoBackFullScreenCloseQuit11.3Metropolis-HastingsMCMC?-:3u?E?d?8?MarkovXn,n=0,1,2,3Xn?G?e?)e?G?Xn+1?!0?Hastings-MetropolisM-H5?E,V4?Markovb?a,?8Ipj,j S,PjSpj=1QG?mS=1,2,N?Markov?=V?Q=(qij)(i,j=1,2,N)MarkovXn,n=0,1,2,?Xn=i)C?PX=j=qij,j=1,2,N25/52kJIkJIGoBackFullScreenCloseQuitXJX=jKXn+1Vij=?jV1 ij=?i.XJ?ij=minjqjiiqij,1KMarkovXn,n=0,1,2,?=V?(pij=qijij,ej 6=ipii=qii+Pk6=iqik(1 ik)iqijij=jqjiji?j 6=i=ipij=jqji?j 6=i ij=jqjiiqijji=1ij=1ji=iqijjqji26/52kJIkJIGoBackFullScreenCloseQuit?Metropolis-HastingM-H 1?E?g(|Xn)?)l?Y?l?=?XlT?J?)?-8I?f.2?Xn+1=(Y,VXn,V1 =(Xn,Y)=minf(Y)g(Xn|Y)f(Xn)g(Y|Xn),1312?)l!U(0,1)?UeU?Xl?N(1,21)Y l?N(2,22),?OY=yX?Vf(x|y)=121p1 2exp(1221(1 2)?x (1+12(y 2)?2)Y=yXE,l?kEX|Y=y=1+12(y 2)V arX|Y=y=21(1 2)nO?X=xY l?32/52kJIkJIGoBackFullScreenCloseQuitGibbs?R“N drop X rho 0.65;mu1 1;mu2 2;sigma1 1;sigma2 s1 sqrt(1 rho2)sigma1;s2 X1,for(iin2:N)+x2 Xi 1,2+m1 mu1+rho sigma1/sigma2 (x2 mu2)+Xi,1 rnorm(1,m1,s1)+x1 Xi,1+m2 mu2+rho sigma2/sigma1 (x1 mu1)+Xi,2 b x dim(x)150002colMeansapplyOx?1=1,2=2?1 colMeans(x)10.96329442.0231725 apply(x,2,mean)10.96329442.023172535/52kJIkJIGoBackFullScreenCloseQuitcovOx?9X?cov(x),1,21,0.9838777 0.51442652,0.51442650.6297255 cor(x),1,21,1.0000000 0.65354762,0.65354761.0000000?1=1,2=0.8,=0.65?1?uy?O?C?36/52kJIkJIGoBackFullScreenCloseQuitplotGibbs?:plot(x,main=00,cex=0.5,xlab=bquote(X1),ylab=bquote(X2),ylim=range(x,2):w?k?5K5A?11-2 Gibbs?:37/52kJIkJIGoBackFullScreenCloseQuit11.5?“dddMCMC?OOOlcA!w?MCMCdnMetropolis 1953cJ)n?pE,OKTL?)-?MarkovTMarkov?-a,?|TMarkov5?1 MCMC?20-V?5?Kcu?ke?r?MCMCk-VO+?)?KL8c MCMC3OA!7K?“dOL?+5?2?A?!?MCMC3?“d?O?A38/52kJIkJIGoBackFullScreenCloseQuit11.5.1?“dddMCMC?OOOA?“d?1K?(J?I?1pOp(JqOm8ck?|?5;?1pOp?MCMC5b?p(y|)?Vy|*=(1,d)?O?p?“d?1?O?k?V()K?V(|y)=p(y|)()Rp(y|)()dp(y|)()(11.5.1)39/52kJIkJIGoBackFullScreenCloseQuit?3SK?11.5.1?E,?/|?$D?Ak5?151?,?(JMCMC5)d?53?cL?e(g+1)?6(g)?=5LK(g),A|y)=P(g+1)A|y,(g)?K(,A|y)?(i+1)K(i),|y)kU48?S?(i+1)?E?8I,K?cv?n0gS“?G(n0+1),(n0+2),(n0+G)5g(|y)?Gv?K?(n0+1),(n0+2),(n0+G)40/52kJIkJIGoBackFullScreenCloseQuit?“nsurrogate bh()=G1GXg=1h(n0+g)(11.5.2)5?Oh()?bl=G1GXg=1(n0+g)l)(11.5.3)5?O?1l.41/52kJIkJIGoBackFullScreenCloseQuitMetropolis-HastingsM-H?gll?qB?=q(,|y)?),?:,?=(?MarkovkX?(?48I pq(,|y)eN?M-H-(g)?ce(g+1)d?)Step1z(0)Step23zgS“g(g=0,1,n0+G)lq(g),|y)?)O(g),|y)=min(1,(|y)(g)|y)q(,(g)|y)q(g),|y)(11.5.4)42/52kJIkJIGoBackFullScreenCloseQuit-(g+1)=(,V(g),|y)(g),V1 (g),|y)Step3DKcn0gS“;?Gg?(n0+1),(n0+2),(n0+G)3SA,|y?eA?1Random walk M-H:?c=+?Xq N(0,V)d(,|y)=min?1,(|y)(|y)?43/52kJIkJIGoBackFullScreenCloseQuit2Independence M-H:-q(,|y)=q(,|y)(,|y)=min?1,(|y)(|y)q(|y)q(|y)?3p?nk7?g?)gb?(1,2)3KK?5-q1(1,1|y,2)q2(2,2|y,1)L?1(1,1|y,2)=min(1,(1|y,2)(1|y,2)q1(1,1|y,2)q1(1,1|y,2)(11.5.5)44/52kJIkJIGoBackFullScreenCloseQuit2(2,2|y,1)=min(1,(2|y,1)(2|y,1)q2(2,2|y,1)q2(2,2|y,1)(11.5.6)?(1|y,2),(1|y,2)?d Bayesn?X(1|y,2)(1,2|y)d3(11.5.5)(11.5.6)?=V?(1,2|y)L?Kz?45/52kJIkJIGoBackFullScreenCloseQuit3 M-H Step1:)1 q1(1(g),1|y,(g)2)V1(g)1,1|y,(g)2)#(g+1)1.Step2:#1?)2 q2(2(g+1),2|y,(g+1)1)V2(g)2,2|y,(g+1)1)#(g+1)2.46/52kJIkJIGoBackFullScreenCloseQuit2Gibbs?3/ezk?/?q1(1(g),1|y,(g)2)=(1|y,(g)2)q2(2(g+1),2|y,(g+1)1)=(2|y,(g+1)1)N?#V1(1(g),1|y,(g)2)=1 A?M-HGibbs?Gibbs?zgS“g(g=1,2,n0+G)l(1|y,(g)2)(g+1)1 l(2|y,(g+1)1)(g+1)2(n0+1),(n0+2),(n0+G)O?ObhG=G1PGg=1h(n0+g).47/52kJIkJIGoBackFullScreenCloseQuit3?3?K?U?IM-HGibbs5?ukO/?uO?M-H48/52kJIkJIGoBackFullScreenCloseQuit11.5.2MCMC?555MCMC?16uMCMC?5e?-(J5gu Tierney(1994).nnn 11.5.1 HHH555b?(g)?kC(|y)K(|y)?CXJ?!?K?8A|P(g)A|y,(0)ZA(|y)|0,a.s.,?G (11.5.7)XJH?KkvZ|h()|(|y)d?h()k49/52kJIkJIGoBackFullScreenCloseQuitbhGZ|h()|(|y)d,a.s.,?G (11.5.8)p?ObhG:=G1PGg=1h(g)n11.5.1J?MCMC?|=?5?!?V9?O5?n?f?k2?A50/52kJIkJIGoBackFullScreenCloseQuitnnn 11.5.2 HHH555%444nnnb?(g)HkC(|y)KkvRh()|2(|y)d?h()kG(bhG Eh)N(0,2h)(11.5.9)Eh=Zh()(|y)d 2h=V ar(h(1)1)+2Xj=2Cov(h(1)1),h(j)1)51/52kJIkJIGoBackFullScreenCloseQuitn11.5.2J?OEh?MCMC?O?ObhG?O?gIOV ar(bhG)=2hG V ar(bhG)?Numerical standard error Monte Carlo standard error(MCSE)LRipley(1987)?1batch means?:-Zg=h(g)(g=1,2,G)rZ1,Z2,ZGyp-U?m?k mBi=m1(Z(i1)m+1+Zim),i=1,2,km?I?y1?gX?u0.05PB:=1kPki=1Bi=bhGB?V ar(B)=1k(k1)Pki=1(Bi B)2?kmGO2hG?O52/52kJIkJIGoBackFullScreenCloseQuit?OV ar(bhG)O-?Ok?f(Inefficiency Factor,Ineff):Ineff(bhG)=V ar(bhG)s2/G=2h/Gs2/Gs2Zg?k?(effective sample size,ESS)5)Ineff?ESS(bhG)=GIneff(bhG)uIneff?,?dL:Ineff(bhG)=Xg=g=1+2Xg=1g(11.5.10)gh()?g?g