70 likes | 519 Views
Metropolis Algorithm Matlab practice. Matlab code taken from Professor Joo -Ho Choi He applied it to with proposal distribution N(x,10). Matlab code for that give in the notes. Here applied to the triangular distribution with U(x-0.25,x+0.25). Metropolis-Hastings algorithm.
E N D
Metropolis Algorithm Matlab practice • Matlab code taken from Professor Joo-Ho Choi • He applied it to with proposal distribution N(x,10). • Matlab code for that give in the notes. • Here applied to the triangular distribution with U(x-0.25,x+0.25)
Metropolis-Hastings algorithm • Metropolis algorithm • A proposal distribution q(x*|x) is symmetric w.r.t x* and x. Then the ratio is simplified. Example: normal pdf at x* with mean at x equals to the vice versa. • Practice with matlab • Generate samples of this distribution using a proposal pdf • As the random walk progresses, the number of samples are increased, and the distribution converges to the target distribution.
Triangular distribution • p=2x in [0,1] • Matlab realization p=@(x) heaviside(x).*heaviside(1-x)*2.*x xx=linspace(-1,2,301); fp=p(xx); plot (xx,fp)
Uniform proposal distribution • q(x)=U(x-0.25,x+0.25) • Sampling: clear; X(1)=0; N=1e4; delta=0.25; for i=1:N-1; x=X(i); xs=x+2*(rand-0.5)*delta u=rand if u<min(1,p(xs)/(p(x)+1.e-10)); X(1+i)=xs; else; X(1+i)=x; end; end • Plotting N0=1; xx=linspace(0,1,26);dx=0.04; %It often pays to have larger N0 nb=histc(X(N0+1:N),xx); bar(xx+dx/2,nb/(N-N0)/dx);
Compare CDFs ecdf(X); hold on xx=linspace(0,1,101); xx2=xx.^2; plot(xx,xx2,'r') legend('ecdf','exact')
practiceproblems • Try triangular distribution with with a normal proposal distribution with different mean and standard deviations. • Do Professor Choi’s example with different starting points. Source: Smithsonian Institution Number: 2004-57325