function [smpl,accept]=mhsample(start,nsamples,varargin) delta=.5; pdf=@normpdf; proppdf=@(x,y)unifpdf(y-x,-delta,delta); proprnd=@(x)x+rand(size(x))*2*delta-delta; sym=false; K=0;%burnin m=1;%thin for k=1:2:length(varargin) if ischar(varargin{k}) switch lower(varargin{k}) case "pdf" if isa(varargin{k+1},"function_handle") pdf=varargin{k+1}; else error("mhsample: pdf must be a function handle"); endif case "proppdf" if isa(varargin{k+1},"function_handle") proppdf=varargin{k+1}; else error("mhsample: proppdf must be a function handle"); endif case "proprnd" if isa(varargin{k+1},"function_handle") proprnd=varargin{k+1}; else error("mhsample: proprnd must be a function handle"); endif case "symmetric" if isa(varargin{k+1},"logical") sym=varargin{k+1}; else error("mhsample: sym must be true or false"); endif case "burnin" if varargin{k+1}>0 K=varargin{k+1}; else error("mhsample: K must be greater than or equal to 0"); endif case "thin" if varargin{k+1}>1 m=varargin{k+1}; else error("mhsample: m must be greater than or equal to 1"); endif case "nchain" warning("mhsample: nchain option not implimented, option ignored"); otherwise warning(["mhsample: Ignoring unknown option " varargin{k}]); endswitch else error("mhsample: Incorect input."); endif endfor ##tic smpl=zeros(nsamples,size(start,2)); smpl(1,:)=start; cx=start; accept=0; i=1; rnd=rand(nsamples*m+K,1); for k=1:nsamples*m+K canacc=rem(k-K,m)==0; px=proprnd(cx); if sym A=(pdf(px))./(pdf(cx)); else A=(pdf(px).*proppdf(cx,px))./(pdf(cx).*proppdf(px,cx)); endif if rnd(k)K accept++; end elseif canacc smpl(i,:)=cx; endif if k>K&&canacc i++; endif end accept./=(nsamples*m); ##toc ##tic ##[smpl acc]=mhsample_ccc(start,nsamples,pdf,proppdf,proprnd,sym,K,m); ##toc end ##%!demo ##%! ## Define function to sample ##%! d=2; ##%! mu=[-1;2]; ##%! Sigma=rand(d); ##%! Sigma=(Sigma+Sigma'); ##%! Sigma+=eye(d)*abs(eigs(Sigma,1,"sa"))*1.1; ##%! pdf=@(x)(2*pi)^(-d/2)*det(Sigma)^-.5*exp(-.5*sum((x.'-mu).*(Sigma\(x.'-mu)),1)); ##%! ##Inputs ##%! start=ones(1,2); ##%! nsamples=500; ##%! sym=true; ##%! K=500; ##%! m=10; ##%! [smpl,accept]=mhsample(start,nsamples,"pdf",pdf,"symmetric",sym,"burnin",K,"thin",m); ##%! figure;hold on; ##%! plot(smpl(:,1),smpl(:,2),'x'); ##%! [x,y]=meshgrid(linspace(-6,4),linspace(-3,7)); ##%! z=reshape(pdf([x(:) y(:)]),size(x)); ##%! mesh(x,y,z,"facecolor","None"); ##%! ## Using sample points to find the volume of half a sphere with radius of .5 ##%! f=@(x)((.25-(x(:,1)+1).^2-(x(:,2)-2).^2).^.5.*(((x(:,1)+1).^2+(x(:,2)-2).^2)<.25)).'; ##%! int=mean(f(smpl)./pdf(smpl)) ##%! errest=std(f(smpl)./pdf(smpl))/nsamples^.5 ##%! trueerr=abs(2/3*pi*.25^(3/2)-int) ##%! mesh(x,y,reshape(f([x(:) y(:)]),size(x)),"facecolor","None"); ##%! hold off; ##%!demo ##%! ##Integrate truncated normal distribution to find normilization constant ##%! pdf=@(x)exp(-.5*x.^2)/(pi^.5*2^.5); ##%! nsamples=1e3; ##%! [smpl,accept]=mhsample(1,nsamples,"pdf",pdf,"symmetric",true,"thin",4); ##%! f=@(x)heaviside(x+2).*exp(-.5*x.^2).*(1-heaviside(x-2)); ##%! int=mean(f(smpl)./pdf(smpl)) ##%! errest=std(f(smpl)./pdf(smpl))/nsamples^.5 ##%! trueerr=abs(erf(2^.5)*2^.5*pi^.5-int)