function [smpl,accept]=mhsample(start,nsamples,varargin) delta=.5; sizestart=size(start); pdf=@(x)mvnpdf(x,zeros(1,sizestart(2)),eye(sizestart(2))); proppdf=@(x,y)prod(unifpdf(y-x,-delta,delta),2); logpdf=[]; logproppdf=[]; proprnd=@(x)x+rand(size(x))*2*delta-delta; sym=false; K=0;%burnin m=1;%thin nchain=1; 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 "logpdf" if isa(varargin{k+1},"function_handle") pdf=varargin{k+1}; else error("mhsample: pdf must be a function handle"); endif case "logproppdf" 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" if varargin{k+1}>=1 nchain=varargin{k+1}; else error("mhsample: nchain must be greater than or equal to 1"); endif otherwise warning(["mhsample: Ignoring unknown option " varargin{k}]); endswitch else error(["mhsample: " varargin{k} " is not a valid property."]); endif endfor ##tic if ~isempty(pdf) && isempty(logpdf) logpdf=@(x)rloge(pdf(x)); end if ~isempty(proppdf) && isempty(logproppdf) logproppdf=@(x,y)rloge(proppdf(x,y)); end if length(sizestart)==2 sizestart=[sizestart 0]; end smpl=zeros(nsamples,sizestart(2),nchain); if all(sizestart([1 3])==[1 nchain]) smpl(1,:,:)=start; elseif all(sizestart([1 3])==[nchain 0]) smpl(1,:,:)=permute(start,[3, 2, 1]); elseif all(sizestart([1 3])==[1 0]) smpl(1,:,:)=repmat(start,[1, 1, nchain]); end cx=permute(smpl(1,:,:),[3 2 1]); accept=zeros(nchain,1); i=1; rnd=log(rand(nchain,nsamples*m+K)); for k=1:nsamples*m+K canacc=rem(k-K,m)==0; px=proprnd(cx); if sym A=(logpdf(px))-(logpdf(cx)); else A=(logpdf(px)+logproppdf(cx,px))-(logpdf(cx)+logproppdf(px,cx)); endif ac=rnd(:,k)K&&canacc i++; endif end accept./=(nsamples*m+K); ##toc ##tic ##[smpl acc]=mhsample_ccc(start,nsamples,pdf,proppdf,proprnd,sym,K,m); ##toc end function y=rloge(x) y=-inf(size(x)); xg0=x>0; y(xg0)=log(x(xg0)); 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)