monte carlo控制和加速收敛方法.ppt

控制和加速收敛（续） ;知识点回顾;h=function(x){(cos(50*x)+sin(20*x))^2} x=matrix(h(runif(200*10^4)),ncol=200) estint=apply(x,2,cumsum)/(1:10^4) plot(estint[,1],ty=l,col=blue,ylim=c(.8,1.2)) y=apply(estint,1,quantile,c(.025,.975)) polygon(c(1:10^4,10^4:1),c(y[1,],rev(y[2,])),col=wheat) boot=matrix(sample(x[,1],200*10^4,rep=T),nrow=10^4,ncol=200) bootit=apply(boot,2,cumsum)/(1:10^4) bootup=apply(bootit,1,quantile,.975) bootdo=apply(bootit,1,quantile,.025) polygon(c(1:10^4,10^4:1),c(bootup,rev(bootdo)),col=gold) 运行结果如下图：;;2.Asymptotic variance of importance sampling estimators （重要性抽样估计的渐进方差）;原积分是：;当 x 是独立同分布（iid）时，它的方差是：;下面我们看一个实际的例子：;运行结果见下图：;3.Effective sample size and perplexity （有效的取样大小和perplexity）;第二种估计取样大小的方法：;用第一种方法的程序： norma=matrix(rnorm(500*10^4),ncol=500)+2.5 weit=1/(1+norma^2) ess=apply(weit,2,cumsum)^2/apply(weit^2,2,cumsum) cocha=matrix(rcauchy(500*10^4),ncol=500) wach=dnorm(cocha,mean=2.5) wachd=wach ech=apply(wachd,2,cumsum)^2/apply(wachd^2,2,cumsum) plot(ess,col=red,xlim=c(0,10000)) plot(ech,col=blue,xlim=c(0,10000) 用第二种方法的程序： sumweit=apply(weit,2,cumsum) plex=(apply(weit*log(weit),2,cumsum)/sumweit)-log(sumweit) chumweit=apply(wachd,2,cumsum) plech=(apply(wachd*log(wachd),2,cumsum)/chumweit)- + log(chumweit) plob=apply(exp(plex),1,quantile,c(.025,.975)) ploch=apply(exp(plech),1,quantile,c(.025,.975));第一种方法的结果如下：;第二种方法的结果;4. Simultaneous monitoring（同步监控）;如果我们考虑以下积分;结果如下图：;5. Rao—Blackwellization and deconditioning;下面我们具体看看这种方法 ;一个具体的例子如下：;用以下程序产生一个（xi,yi）的样本：y=sqrt(rchisq(Nsim,df=nu)/nu)x=rnorm(Nsim,mu,sigma/y);计算程序为：d1=cumsum(exp(-x^2))/(1:Nsim)d2=cumsum(exp(3^2/(1+2*(0.5/y)^2))/sqrt(1+2*(0.5/y)^2))/(1:Nsim);当使用一个混合分布作为重要函数，;我们再来看一种实现Rao-Blackwellization的特殊形式;;6.Acceleration methods（加速方法）;6.1.Correlated simulations（相关模拟）;;;;程序如下：; 在更一般的情况下，在两个风险函数之间创建一个足够强的相关性是比较困难的，当???个函数之间存在简单的转换关系，如尺度转换，我们就可以应用统一的样本，否则就会造成一些全局的负面影响。比如说，如果两个函数都是通过取舍法来模拟取样，那么因为获得的随机变量的个数是随机的，这时这两个函数就不能应用相同的样本。;6.2.Antithetic variables（对偶变量）;uref=runif(10^4) x=h(uref)