实例介绍
【实例简介】
是R中mcmc算法的一个例子和介绍，对新接触R与mcmc算法的人有很重要的参考价值
[1]0.008 The arguments to the metrop function here (there are more we don't use here) are an R function(here lupost that evaluates the log unnormalized density of the desired stationary distribution (here a posterior distribution) of the Markov chain. Note that( although this example does not exhibit the phenomenon) that the unnormalized density may be zero, in which case the log unnormalized density is -Inf an initial state(here beta init) of the Markov chain a number of batches(here 1e3) for the Markov chain. This combines with batch length and spacing(both 1 by default) to determine the number of iterations done additional arguments(here x and y) supplied to provided functions(here luposa there is no"burn-in"argument. although burn-in is easily accomplished if desired(morc on this below The output is in the component outsbatch returned by the metrop function Wo'll look at it presently, but first we nccd to adjust the proposal to gct a highcr acceptance rate(out saccept). It is generally accepted( Gelman, Roberts, and Gilks, 1996)t, hat an acceptance rate of about 20% is right, although this rec- ommendation is based on the asymptotic analysis of a toy problem (simulating a multivariate normal distribution) for which one would never use MCMC and is very unrepresentative of difficult MCMC applications Geyer and Thompson(1995) came to a similar conclusion, that a 20% accep tance ratc is about right, in a vcry diffcrent situation. But thcy also warned that d 20% acceptance rate could be very wrong and produced all exanple where i 20% acceptance rate was impossible and attempting to reduce the acceptance rate below 70% would keep the sampler from ever visiting part of the state space So the 20%o magic number must be considered like other rules of thumb we teach in intro courses (like n> 30 meansmeans normal approximation is valid). We know these rules of thumb can fail. There are examples in the literature where they do fail. We keep repeating them because we want something simple to tell beginners, and they are all right for some problems Be that as it may we try for 20% out <-metrop(out, scale =0 y =y utsaccept [1]0.739 t <-metrop(out, scale=0 outaccept 3 [1]0.371 out <-metrop (out, scale =0.5, x=x, y= y) > outaccept [1]0.148 out < metro(out, scale =0.4, x outsaccept [1]0.209 Here the first, argument to each instance of the metro function is the output f a previous invocation. The Markov chain continues where the previous run stopped doing just what it would have done if it had kept going the initial state and random seed being the final state and random seed of the previous invocation. Everything stays the same except for the arguments(here scale that are supplied) The argument scale controls the size of the Metropolis "normal random walk" proposal. The default is scale =1. Big steps give lower acceptance rates. Small steps give higher. We want something about 20%. It is also possible to make scale a vector or a matrix. See help(metro) Because each run starts where the last one stopped (when the first argument to metro is the output of the previous invocation), each run serves as"burn-in for its successor(assuming that any part of that run was worth anything at all) 3 Diagnostics O.K. That does it for the acceptance rate. So let s do a longer run and look .t the result out <-metrop (out, batch =10000,x=x, y=y) outsaccept [1]0.2345 > outstime user system elapsed 0.4000.0000.403 Figure 1 (page 5) shows the time series plot made by the R statement plot(ts (out batch)) Another way to look at the output is all autocorrelation plot. Figure 2 (page 6)shows the time series plot made by the R statement Figure 1: Timle series plot of MCMC output. ud叫 TTTTTT ATTTIT TT Figure 2: AutocorrelatiOn plot of MCMC output 6 >acf (outmatch) As with any multiplot plot, these are a bit hard to read. Readers are invited to make the scparatc plots to gct a bctter picturc. As with all" diagnostic" plots in MCMC, these dont"diagnose subtle problems. As http://www.stat.umn.edu/charlie/mcmc/diag.htm1 says The purpose of regression diagnostics is to find obvious, gross embarrassing problems that jump out of simple plots The time series plots will show obvious nonstationarity. They will not show nonobuious nonstationarity. Thcy provide no guarantee whatsoever that your Markov chain is sanpling anything remotely resembling the correct stationary distribution(with log unnormalized density lupost). In this very easy problem we do not expect any convergence difficulties and so believe what the diagnostics seem to show, but one is a fool to trust such diagnostics in difficult problems The autocorrelation plots seem to show that the the autocorrelations are negligible after about lag 25. This diagnostic inference is reliable if the sampler is actually working(has nearly reached equilibrium)and worthless otherwise Thus batches of length 25 should be sufficient, but let's use length 100 to be 4 Monte carlo estimates and standard errors out <-metrop(out, batch =100, blen =100, outfun function(z, .)c(z, z2),x=x, y=y, outsaccept [1]0.2332 outstime user system elapsed 0.4400.0000.437 We have added an argument outfun that gives the"functional "of the state we want to average. For this problem we are interested in both posterior mean and variance. Mean is easy, just average the variables in question. But variance is a little tricky. We nced to usc the identity (X)=E(X2)-E(X to write variance as a function of two things that can be estimated by simple averages. Hence we want to average the state itself and the squares of each component. Hence our outfun returns c(z, z"2) for an argument (the state tor )z The... arguinent to outfun is required, since the functioN is also pass the other arguments(here x and y)to metro. mple means The grand means(means of batch means)are apply (out batch, 2, mean) [1]0.65319500.79203421.17010750.50773310.7488265 [6]0.51457510.75607751.49738070.39138370.7244162 The first 5 numbers are the Monte Carlo estimates of the posterior means. The second 5 numbers are the Monte Carlo estimates of the posterior absolute second moments. We get the posterior variances by apply (outmatch, 2, >mu <-foo[1: 5] sigmasq < foo[6: 10]-mu"2 mu [1]0.65319500.79203421.17010750.50773310.7488265 [1]0.087911340.128759240.128229240.133590810.16367507 Monte Carlo st andard errors(Cse) are calculated from the batch means This is simplest for the means mu mcse < apply(out batch, 1: 51, 2, sd)/sqrt(out Sabat ch) mu. mcse [1]0.012242600.014179160.017931290.014685940.01582040 The extra factor sart(outsnbatch) arises because the batch means have vari- ance o2/b where b is the batch length, which is outsblen, whereas the overall mcans mu havc variance a2/n whcrc n is thc total numbcr of itcrations, which outblen outonbatch 4.2 Functions of means To get the mose for the posterior variances we apply the delta method. Let ui denote the sequence of batch means of the first kind for one parameter and u the grand mean(the estimate of the posterior mean of that parameter). let v, denotc the scqucncc of batch mcans of the sccond kind for the samc parameter anld v the grand meall(the estimate of the posteric eult of that parameter), and let u= E(a and 1=E(u). Then the delta method linearizes the nonlinear function 9(u, v) as △9(1,D) saying that nas the same asymptotic normal distribution as which of course, has variance 1/ batch times that of (v2-y)-21(-1) and this variance is estimated by ∑[(v-7)-2a(u;-a) batch DO u <-outsbatch[, 1:51 v < outsbatchL, 6: 10] ubar <-apply (u, 2, mean) bar <-apply(v, 2, mean) delta < sweep(u, 2, ubar) deltav < sweep foo <-sweep(deltan, 2, ubar, *" sigmasq mcse < sqrt (apply((deltav -2* foo)"2 2, mean/out snatch) sigmasg mcse [1]0.0042416370.0072922240.0072713900.007374833 [5]0.008175832 does the MCse for the posterior variance Let's just check that this complicated sweep and apply stuff does do the right thing sqrt(mean(((vL, 2]-vbarL21)-2 ubar[2]* ul, 2]-ubar[2]))2)/outonbatch [1]0.007292224 Comment Through version 0.5 of this vignette it contained an incorrect pro- cedure for calculating this MCSe, justified by a handwave. Essentially, it said to use the standard deviation of the batch means called v here, which appears to be very conservative 4.3 Functions of Functions of means If we are also interested in the posterior standard deviation(a natural ques tion, although not asked on the exam problem), the delta method gives its standard error in terms of that for the variance sigma < sqrt(sigmasq sigma mcse < sigmasq. mcse/(2* sigma) Igma [1]0.29649850.35883040.35809110.36550080.4045678 > sigma a. mcse [1]0.0071528820.0101611020.0101529890.010088669 [5]0.010104403 5A Final rur So thats it. The only thing left to do is a little more precision(the exam problem directed"use a long enough run of your Markov chain sampler so that the mcse are less than 0.01) out <-metrop(out, batch 500, blen =400, X =X, y y outaccept [1]0.235155 0ut$七ime user system elapsed 8.8240.0008.826 foo < apply (out batch, 2, mean) mu < foo[1: 5] sigmasg < foo 16: 10]-mu"2 mu [1]0.66246500.79410131.17127100.50663260.7261414 sIgmasq [1]0.091892460.133230540.132308110.128712930.15978638 mu mcse < apply (out batch, 1: 51, 2, sd)/sqrt(outsnbatch) mu. mcse [1]0.0029601280.0036474200.0037878550.003632080 [5]0.004273624 【实例截图】
【核心代码】
是R中mcmc算法的一个例子和介绍，对新接触R与mcmc算法的人有很重要的参考价值
[1]0.008 The arguments to the metrop function here (there are more we don't use here) are an R function(here lupost that evaluates the log unnormalized density of the desired stationary distribution (here a posterior distribution) of the Markov chain. Note that( although this example does not exhibit the phenomenon) that the unnormalized density may be zero, in which case the log unnormalized density is -Inf an initial state(here beta init) of the Markov chain a number of batches(here 1e3) for the Markov chain. This combines with batch length and spacing(both 1 by default) to determine the number of iterations done additional arguments(here x and y) supplied to provided functions(here luposa there is no"burn-in"argument. although burn-in is easily accomplished if desired(morc on this below The output is in the component outsbatch returned by the metrop function Wo'll look at it presently, but first we nccd to adjust the proposal to gct a highcr acceptance rate(out saccept). It is generally accepted( Gelman, Roberts, and Gilks, 1996)t, hat an acceptance rate of about 20% is right, although this rec- ommendation is based on the asymptotic analysis of a toy problem (simulating a multivariate normal distribution) for which one would never use MCMC and is very unrepresentative of difficult MCMC applications Geyer and Thompson(1995) came to a similar conclusion, that a 20% accep tance ratc is about right, in a vcry diffcrent situation. But thcy also warned that d 20% acceptance rate could be very wrong and produced all exanple where i 20% acceptance rate was impossible and attempting to reduce the acceptance rate below 70% would keep the sampler from ever visiting part of the state space So the 20%o magic number must be considered like other rules of thumb we teach in intro courses (like n> 30 meansmeans normal approximation is valid). We know these rules of thumb can fail. There are examples in the literature where they do fail. We keep repeating them because we want something simple to tell beginners, and they are all right for some problems Be that as it may we try for 20% out <-metrop(out, scale =0 y =y utsaccept [1]0.739 t <-metrop(out, scale=0 outaccept 3 [1]0.371 out <-metrop (out, scale =0.5, x=x, y= y) > outaccept [1]0.148 out < metro(out, scale =0.4, x outsaccept [1]0.209 Here the first, argument to each instance of the metro function is the output f a previous invocation. The Markov chain continues where the previous run stopped doing just what it would have done if it had kept going the initial state and random seed being the final state and random seed of the previous invocation. Everything stays the same except for the arguments(here scale that are supplied) The argument scale controls the size of the Metropolis "normal random walk" proposal. The default is scale =1. Big steps give lower acceptance rates. Small steps give higher. We want something about 20%. It is also possible to make scale a vector or a matrix. See help(metro) Because each run starts where the last one stopped (when the first argument to metro is the output of the previous invocation), each run serves as"burn-in for its successor(assuming that any part of that run was worth anything at all) 3 Diagnostics O.K. That does it for the acceptance rate. So let s do a longer run and look .t the result out <-metrop (out, batch =10000,x=x, y=y) outsaccept [1]0.2345 > outstime user system elapsed 0.4000.0000.403 Figure 1 (page 5) shows the time series plot made by the R statement plot(ts (out batch)) Another way to look at the output is all autocorrelation plot. Figure 2 (page 6)shows the time series plot made by the R statement Figure 1: Timle series plot of MCMC output. ud叫 TTTTTT ATTTIT TT Figure 2: AutocorrelatiOn plot of MCMC output 6 >acf (outmatch) As with any multiplot plot, these are a bit hard to read. Readers are invited to make the scparatc plots to gct a bctter picturc. As with all" diagnostic" plots in MCMC, these dont"diagnose subtle problems. As http://www.stat.umn.edu/charlie/mcmc/diag.htm1 says The purpose of regression diagnostics is to find obvious, gross embarrassing problems that jump out of simple plots The time series plots will show obvious nonstationarity. They will not show nonobuious nonstationarity. Thcy provide no guarantee whatsoever that your Markov chain is sanpling anything remotely resembling the correct stationary distribution(with log unnormalized density lupost). In this very easy problem we do not expect any convergence difficulties and so believe what the diagnostics seem to show, but one is a fool to trust such diagnostics in difficult problems The autocorrelation plots seem to show that the the autocorrelations are negligible after about lag 25. This diagnostic inference is reliable if the sampler is actually working(has nearly reached equilibrium)and worthless otherwise Thus batches of length 25 should be sufficient, but let's use length 100 to be 4 Monte carlo estimates and standard errors out <-metrop(out, batch =100, blen =100, outfun function(z, .)c(z, z2),x=x, y=y, outsaccept [1]0.2332 outstime user system elapsed 0.4400.0000.437 We have added an argument outfun that gives the"functional "of the state we want to average. For this problem we are interested in both posterior mean and variance. Mean is easy, just average the variables in question. But variance is a little tricky. We nced to usc the identity (X)=E(X2)-E(X to write variance as a function of two things that can be estimated by simple averages. Hence we want to average the state itself and the squares of each component. Hence our outfun returns c(z, z"2) for an argument (the state tor )z The... arguinent to outfun is required, since the functioN is also pass the other arguments(here x and y)to metro. mple means The grand means(means of batch means)are apply (out batch, 2, mean) [1]0.65319500.79203421.17010750.50773310.7488265 [6]0.51457510.75607751.49738070.39138370.7244162 The first 5 numbers are the Monte Carlo estimates of the posterior means. The second 5 numbers are the Monte Carlo estimates of the posterior absolute second moments. We get the posterior variances by apply (outmatch, 2, >mu <-foo[1: 5] sigmasq < foo[6: 10]-mu"2 mu [1]0.65319500.79203421.17010750.50773310.7488265 [1]0.087911340.128759240.128229240.133590810.16367507 Monte Carlo st andard errors(Cse) are calculated from the batch means This is simplest for the means mu mcse < apply(out batch, 1: 51, 2, sd)/sqrt(out Sabat ch) mu. mcse [1]0.012242600.014179160.017931290.014685940.01582040 The extra factor sart(outsnbatch) arises because the batch means have vari- ance o2/b where b is the batch length, which is outsblen, whereas the overall mcans mu havc variance a2/n whcrc n is thc total numbcr of itcrations, which outblen outonbatch 4.2 Functions of means To get the mose for the posterior variances we apply the delta method. Let ui denote the sequence of batch means of the first kind for one parameter and u the grand mean(the estimate of the posterior mean of that parameter). let v, denotc the scqucncc of batch mcans of the sccond kind for the samc parameter anld v the grand meall(the estimate of the posteric eult of that parameter), and let u= E(a and 1=E(u). Then the delta method linearizes the nonlinear function 9(u, v) as △9(1,D) saying that nas the same asymptotic normal distribution as which of course, has variance 1/ batch times that of (v2-y)-21(-1) and this variance is estimated by ∑[(v-7)-2a(u;-a) batch DO u <-outsbatch[, 1:51 v < outsbatchL, 6: 10] ubar <-apply (u, 2, mean) bar <-apply(v, 2, mean) delta < sweep(u, 2, ubar) deltav < sweep foo <-sweep(deltan, 2, ubar, *" sigmasq mcse < sqrt (apply((deltav -2* foo)"2 2, mean/out snatch) sigmasg mcse [1]0.0042416370.0072922240.0072713900.007374833 [5]0.008175832 does the MCse for the posterior variance Let's just check that this complicated sweep and apply stuff does do the right thing sqrt(mean(((vL, 2]-vbarL21)-2 ubar[2]* ul, 2]-ubar[2]))2)/outonbatch [1]0.007292224 Comment Through version 0.5 of this vignette it contained an incorrect pro- cedure for calculating this MCSe, justified by a handwave. Essentially, it said to use the standard deviation of the batch means called v here, which appears to be very conservative 4.3 Functions of Functions of means If we are also interested in the posterior standard deviation(a natural ques tion, although not asked on the exam problem), the delta method gives its standard error in terms of that for the variance sigma < sqrt(sigmasq sigma mcse < sigmasq. mcse/(2* sigma) Igma [1]0.29649850.35883040.35809110.36550080.4045678 > sigma a. mcse [1]0.0071528820.0101611020.0101529890.010088669 [5]0.010104403 5A Final rur So thats it. The only thing left to do is a little more precision(the exam problem directed"use a long enough run of your Markov chain sampler so that the mcse are less than 0.01) out <-metrop(out, batch 500, blen =400, X =X, y y outaccept [1]0.235155 0ut$七ime user system elapsed 8.8240.0008.826 foo < apply (out batch, 2, mean) mu < foo[1: 5] sigmasg < foo 16: 10]-mu"2 mu [1]0.66246500.79410131.17127100.50663260.7261414 sIgmasq [1]0.091892460.133230540.132308110.128712930.15978638 mu mcse < apply (out batch, 1: 51, 2, sd)/sqrt(outsnbatch) mu. mcse [1]0.0029601280.0036474200.0037878550.003632080 [5]0.004273624 【实例截图】
【核心代码】
标签：
好例子网口号：伸出你的我的手 — 分享！
小贴士
感谢您为本站写下的评论，您的评论对其它用户来说具有重要的参考价值，所以请认真填写。
- 类似“顶”、“沙发”之类没有营养的文字，对勤劳贡献的楼主来说是令人沮丧的反馈信息。
- 相信您也不想看到一排文字/表情墙，所以请不要反馈意义不大的重复字符，也请尽量不要纯表情的回复。
- 提问之前请再仔细看一遍楼主的说明，或许是您遗漏了。
- 请勿到处挖坑绊人、招贴广告。既占空间让人厌烦，又没人会搭理，于人于己都无利。
关于好例子网
本站旨在为广大IT学习爱好者提供一个非营利性互相学习交流分享平台。本站所有资源都可以被免费获取学习研究。本站资源来自网友分享，对搜索内容的合法性不具有预见性、识别性、控制性，仅供学习研究，请务必在下载后24小时内给予删除，不得用于其他任何用途，否则后果自负。基于互联网的特殊性，平台无法对用户传输的作品、信息、内容的权属或合法性、安全性、合规性、真实性、科学性、完整权、有效性等进行实质审查；无论平台是否已进行审查，用户均应自行承担因其传输的作品、信息、内容而可能或已经产生的侵权或权属纠纷等法律责任。本站所有资源不代表本站的观点或立场，基于网友分享，根据中国法律《信息网络传播权保护条例》第二十二与二十三条之规定，若资源存在侵权或相关问题请联系本站客服人员，点此联系我们。关于更多版权及免责申明参见 版权及免责申明
网友评论
我要评论