Adapting the sampler\u00a0throughout the chain has the\u00a0potential to be more efficient than fixing it at one value. Imagine sampling values from the conditional distribution of beta when the posterior resembles the plot below.\u00a0Ideally we would like\u00a0to use a proposal distribution for beta with a small standard deviation\u00a0when alpha is around 10 but a proposal distribution with a much larger standard deviation when alpha is around 20. Since we are unlikely to have the information necessary for designing such a sampler when the chain starts, wouldn’t it be helpful if the chain could learn how to adapt the proposal distribution as it goes along?<\/p>\n
A first sight adaptation looks attractive but before we start considering it, there is one very important message:<\/p>\n adapting the sampler during an MCMC run is dangerous.<\/span><\/p>\n So if you are not sure that it is valid, then don’t do it.<\/p>\n Beware the mathematician<\/strong><\/p>\n In order for an MCMC algorithm to converge, we want the chain to have detailed balance, that is, it needs the property that for any two values of the parameter \u03b8, say \u03b81 and \u03b82<\/p>\n P(\u03b81) P(move \u03b81 to \u03b82) = P(\u03b82) P(move \u03b82 to \u03b81)<\/p>\n Typically the moves are made dependent on (i) a proposal distribution and (ii) an acceptance probability, so the requirement for detailed balance becomes,<\/p>\n P(\u03b81) P(propose \u03b81 to \u03b82) P(accept \u03b81 to \u03b82) = P(\u03b82) P(propose \u03b82 to \u03b81) P(accept \u03b82 to \u03b81)<\/p>\n If we adapt the proposal distribution then we alter P(propose \u03b81 to \u03b82) as the chain progresses and we will have great difficulty in choosing an acceptance probability that ensures balance.<\/p>\n A simple solution to the problem is to insist that the changes to the proposal distribution get smaller over time; eventually they will fade away and the proposal distribution will almost be fixed and as a consequence we will eventually be able to choose acceptance probabilities that guarantee detailed balance. This is usually called diminishing adaptation.<\/p>\n Now there are mathematical proofs to the effect that if the changes in the proposal distribution fade towards zero then you can easily make the chain converge to the correct distribution. Unfortunately these are mathematical proofs about what will happen if the chain is run for ever and of course we do not run chains of infinite length.<\/p>\n Instead of accepting the mathematical proof, we must ask ourselves a more practical question. Suppose I run an adaptive chain for 10,000 iterations and maybe the proposal distribution is changed by noticeable, but diminishing amounts; will this have a meaningful impact on my estimated posterior distributions?<\/p>\n Remember that if we restrict adaptation to the discarded burn-in, as I do with the adapt option, then we will not have a problem.<\/p>\n A program for investigating adaptation<\/strong><\/p>\n Here is a relatively general program for a random walk Metropolis-Hasting sampler\u00a0that simulates values from a normal prior distribution with mean 1 and standard deviation 2. This may look rather unusual since we are used to simulating from posterior distributions that depend on data,\u00a0while here we are simulating from a prior.\u00a0It is also true that\u00a0sampling from N(1,sigma=2) is trivial and does not need MCMC but\u00a0there are advantages in testing an algorithm on a problem\u00a0for which\u00a0we know the correct answer.<\/p>\n To try to be clear, I will use sigma to represent the standard deviation of the\u00a0target distribution\u00a0and sd for the standard deviation of the proposal distribution.<\/p>\n I’ve highlighted 3 lines by labelling them as A, B and C. By changing those lines we can alter then properties of the algorithm.<\/p>\n The random walk updates have a standard deviation that is controlled by the local sd<\/strong> and the proposal distribution is altered by changing sd\u00a0by an amount lambda<\/strong>. As the program stands, I set sd=1 and lambda=1 so the proposal distribution does not change. The mean posterior estimate of theta is 1.089 and its standard error (sem) is 0.0939. The standard deviation of the chain (estimating sigma) is 2.006. Here is a trace plot of the chain.<\/p>\n The chain mixes quite slowly because the proposal standard deviation is fixed at a value well below its optimal value. Setting line (A) to sd =4 produces better mixing with an estimate of 1.030 and a sem of 0.0422.<\/p>\n Both of these are valid chains that have detailed balance because sd is fixed. So now let\u2019s allow sd to change and see what happens. If we return line (A) to read sd = 1 but put line (B) to lambda = 1.02. The proposal distribution will broaden when moves are accepted and narrow when a proposal is rejected. Eventually the chain will settle to reject about 50% of proposals. This is the algorithm used during the burn-in by –mcmcrun<\/strong>– whenever the adapt option is set.<\/p>\n To create a diminishing effect we can set line (C) to be lambda = 1 + 0.02*0.999^`i\u2019. Now we can see the impact of the diminishing changes in the way that sd settles around 4 once we pass iteration 4000 and lambda gets close to 1. The algorithm has the flexibility to adjust but we soon learn about the sd and the changes to become negligible. The results are fine with a posterior mean is 0.973 and sem=0.0414.<\/p>\n![\"contour\"](\"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/files\/2014\/11\/contour.png\")
<\/a><\/p>\n clear<\/span>\r\n set seed 51817<\/span>\r\n set obs 10000<\/span>\r\n *----------------------------<\/span>\r\n * INITIAL VALUES<\/span>\r\n *----------------------------<\/span>\r\n \/*A*\/ local sd = 1<\/span>\r\n \/*B*\/ local lambda = 1<\/span>\r\n gen sd = `sd'<\/span>\r\n gen lambda = `lambda'<\/span>\r\n gen theta = 1<\/span>\r\n *----------------------------<\/span>\r\n * 10,000 iterations<\/span>\r\n *----------------------------<\/span>\r\n local oldlogL = log(normalden(1,1,2))<\/span>\r\n forvalues i=2\/10000 {<\/span>\r\n qui replace sd = `sd' in `i'<\/span>\r\n qui replace lambda = `lambda' in `i'<\/span>\r\n local j = `i' - 1<\/span>\r\n *-------------------------------------<\/span>\r\n * METROPOLIS-HASTINGS<\/span>\r\n * & UPDATE PROPOSAL DISTRIBUTION (SD)<\/span>\r\n *-------------------------------------<\/span>\r\n local newtheta = theta[`j']+rnormal(0,`sd')<\/span>\r\n local logL = log(normalden(`newtheta',1,2))<\/span>\r\n if log(runiform()) < `logL' - `oldlogL' {<\/span>\r\n qui replace theta = `newtheta' in `i'<\/span>\r\n local oldlogL = `logL'<\/span>\r\n local sd = `lambda'*`sd'<\/span>\r\n }<\/span>\r\n else {<\/span>\r\n local sd = `sd'\/`lambda'<\/span>\r\n qui replace theta = theta[`j'] in `i'<\/span>\r\n }<\/span>\r\n *----------------------------<\/span>\r\n * UPDATE DIMINISHING FACTOR<\/span>\r\n *----------------------------<\/span>\r\n \/*C*\/ local lambda = `lambda'<\/span>\r\n qui replace lambda = `lambda' in `i'<\/span>\r\n }<\/span>\r\n histogram theta , \/\/\/<\/span>\r\n addplot( function y = normalden(x,1,2), range(-5 7) ) \/\/\/<\/span>\r\n leg(off)<\/span>\r\n mcmcstats theta<\/span>\r\n count if theta != theta[_n-1]<\/span>\r\n mcmctrace theta lambda sd , cgopt(row(3))<\/span><\/strong><\/pre>\n<\/a><\/p>\n<\/a>
\nThe proposal standard deviation moves quickly towards the (3,5) range. The mean of theta is 1.035 with a sem of 0.0416. The changes in the proposal distribution have no material impact the convergence. Sometimes sd is too large and sometimes it is too small but the effect is small and seems to cancel.<\/p>\n