Metropolis hasting algorithm

Metropolis hasting algorithm:

The MH algorithm is also a simulation-based MCMC technique, and a very important application is sampling from a given probability distribution. The main principle is to construct a subtle Markov chain, so that the steady state of the chain is your given probability density. Its merits, needless to say, are naturally capable of dealing with the complex probability density of mathematical forms. Some say that single-dimension's MH algorithm with Gibbs sampler almost "invincible".

Today's test process found that the MH algorithm to use the good is also not simple, the transfer parameters inside the set is not very good to get. Even with the simplest Gaussian drift term, the determination of variance is also a headache, requiring different problems to treat, more trials several times. Of course you can always choose the "ideal" number of references.

Or the last mixed Gaussian distribution to do the simulation, the number of simulations is 500,000 times, the probability of the distribution approximation of the degree of such as. Although several obvious "peaks" have come out, there is a very big difference in numbers. It is expected that my drift variance is not selected well. Feel or inverse sampling use, the number of iterations not very much, you can reach a considerable degree of approximation.

Tried the MH algorithm,

R Code:

P=function (X,U1,SIG1,U2,SIG2) {
(1/3) * (1/(sqrt (2*PI) *15) *exp ( -0.5* (x-70) ^2/15^2) +1/(sqrt (2*PI) *11) *exp ( -0.5* (x+80) ^2/11^2) +1/(sqrt (2*PI) *sig1) *exp ( -0.5* (X-U1) ^2/sig1^2) +1/(sqrt (2*PI) *sig2) *exp ( -0.5* (X-U2) ^2/sig2^2))
}


Mh=function (x0,n) {
X=null
X[1] = x0
For (i in 1:n) {
x_can= X[i]+rnorm (1,0,3.25)
D= p (x_can,10,30,-10,10)/P (x[i],10,30,-10,10)
Alpha= min (1,d)
U=runif (1,0,1)
if (U<alpha) {
X[i+1]=x_can}
else{
X[i+1]=x[i]
}
if (round (i/100) ==i/100) print (i)
}
X
}
Z=MH (10,99999)
Z=Z[-10000]
A=seq ( -100,100,0.2)

Plot (density (z), col=1,main= ' estimated density ', ylim=c (0,0.02), Lty=1)
Points (A, p (a,10,30,-10,10), pch= '. ', col=2,lty=2)
Legend (60,0.02,c ("True", "Sim (MH)"), Col=c ((), lty=c)

Metropolis hasting algorithm