tag:blogger.com,1999:blog-7931813.post4711985487589400360..comments2013-06-26T19:34:08.566-04:00Comments on Random Fluctuations: Bayesian analysis: Comparing algorithms Part 1?Avihttp://www.blogger.com/profile/17021838649882950034noreply@blogger.comBlogger5125tag:blogger.com,1999:blog-7931813.post-72746894937615126092011-10-06T21:52:16.749-04:002011-10-06T21:52:16.749-04:00If you're trying to learn Bayesian data analys...If you're trying to learn Bayesian data analysis using MCMC methods, I can happily recommend a book. (I happen to be the author of the book :-) but apparently other people seem to like it too.) Here are links to the book's blog and web page:<br />http://doingbayesiandataanalysis.blogspot.com/<br />http://www.indiana.edu/%7Ekruschke/DoingBayesianDataAnalysis/<br />Happy Bayesian analyzing!John K. Kruschkehttps://www.blogger.com/profile/17323153789716653784noreply@blogger.comtag:blogger.com,1999:blog-7931813.post-15169310440969325842011-08-24T14:09:52.386-04:002011-08-24T14:09:52.386-04:00I usually use the metrop() function I mentioned ab...I usually use the metrop() function I mentioned above, but for some special applications I've written my own samplers (e.g., to sample different groups of variables in different blocks, each with their own proposal distribution). I've also worked with more advanced Monte Carlo algorithms (like Hamiltonian Monte Carlo, which uses gradient information to speed up convergence ... assuming that it's easy and fast to compute partial derivatives of your posterior).<br /><br />By the way, for future reference, the chains you plotted for the Metropolis algorithm are symptomatic of either a too-small step size or a correlated posterior or both.<br /><br />What you see in the traces is that the algorithm is taking very small steps to gradually approach a good region of probability space. That either means the step size is too small and can be simply increased, or else the posterior distribution is so correlated that with independent steps you <i>have</i> to take small steps.<br /><br />The two solutions I mentioned are the typical remedies: either increase the step size, or take correlated steps. If there are a lot of correlated variables some people like to break them into small blocks of related parameters (as I mentioned above) and try to move them separately, so a failure to move in one block won't affect the parameters in other blocks.<br /><br />By contrast, the graphical diagnostic of a step size that's too <em>large</em> is a chain that looks very "stepwise": you get lots of rejections, so the chain stays in the same place for a long time before finally moving.<br /><br />One rule of thumb is to adjust the step sizes until you get your acceptance rate to about 20% or so. (This isn't a hard rule, and the effective sample size is a better measure of how much information is in your chain.) If you have a small acceptance rate, take smaller steps; if you have a high acceptance rate, take bigger steps. In your original example, I think it was producing around a 70% acceptance rate, suggesting that you could take bigger steps in at least one of the variables.Nathan Urbanhttp://www.princeton.edu/~nurban/noreply@blogger.comtag:blogger.com,1999:blog-7931813.post-30587226378066009772011-08-24T12:45:47.741-04:002011-08-24T12:45:47.741-04:00Thank you very much, Dr. Urban. I will definitely ...Thank you very much, Dr. Urban. I will definitely try your suggestions! I agree, the existing packages for R will work much better than my home-brewed code, but I do feel that I have a somewhat better understanding of the how the algorithms are supposed to work now that I have made some attempt at writing them from scratch.<br /><br />In your research, do you use the existing packages, or have you written custom-tailored code for your use?<br /><br />Thank you again!Avihttps://www.blogger.com/profile/17021838649882950034noreply@blogger.comtag:blogger.com,1999:blog-7931813.post-4283995663483928802011-08-24T11:48:50.064-04:002011-08-24T11:48:50.064-04:00Actually, the simplest thing you could do to impro...Actually, the simplest thing you could do to improve performance is just to make the theta step size much bigger. Try sd=10000 instead of sd=100. You can see from the Gibbs plot that the spread of theta is of order 10^4. (Again, in an adaptive sense, you could compute the standard deviation of the thetas you got in the first poorly converged Metropolis run, and use that as the standard deviation of your proposal distribution.)Nathan Urbanhttp://www.princeton.edu/~nurban/noreply@blogger.comtag:blogger.com,1999:blog-7931813.post-9983030954277252792011-08-24T11:39:15.714-04:002011-08-24T11:39:15.714-04:00The Metropolis algorithm isn't performing well...The Metropolis algorithm isn't performing well here because the joint posterior for alpha and theta is highly correlated (plot alpha against theta to see this). Your proposal distribution is proposing moves for alpha and theta that are independent of each other, but it would be better to propose moves in a correlated direction along the high probability surface.<br /><br />In this case, it's probably better to just do Gibbs, but in cases where you can't write down the necessary exact conditional distributions, you can try an adaptive Metropolis approach. Compute the covariance of the poorly converged chain and use it as the proposal distribution to compute second, better converged chain. For example, once you've computed alpha and theta as above, you can compute the covariance of this "preliminary" chain:<br /><br />precov = cov(cbind(alpha,theta))<br /><br />Then run the Metropolis algorithm again. Except this time, in your proposal distribution, propose a new alpha and a new theta at the same time, drawn from a (correlated) bivariate normal distribution centered on the existing point:<br /><br />new.point = rmvnorm(1, mean=c(alpha[i],theta[i]), Sigma=precov)<br />alpha[i+1] = new.point[1]; theta[i+1] = new.point[2]<br /><br />(The 'rmvnorm' function may be found in the 'mvtnorm' package, or you can implement multivariate normal sampling yourself.)<br /><br />It's not always optimal to have the proposal distribution be identical to the covariance of the target posterior, so you may have to scale the 'precov' matrix (e.g., multiply it by 0.5) to get a good acceptance rate (around 25%).<br /><br />This won't work perfectly, since a multivariate normal proposal distribution doesn't look like the actual posterior (which is curved in alpha-theta space), but it will work much better than independent proposals. (I tried it on your problem and it looks a lot closer to the Gibbs output.)<br /><br />See Roberts and Rosenthal, "Examples of adaptive MCMC" for more advanced methods.<br /><br />Also, I know this was a teaching example for yourself, but you can get faster performance by using an R package for MCMC like the 'metrop' function in the 'mcmc' package.Nathan Urbanhttp://www.princeton.edu/~nurban/noreply@blogger.com