From: Brendan on
I can run the example case of generating from a normal distribution fine, and I am now trying to generate from an arbitrary multivariate distribution, which does not seem to be working. Each variable seems to be at the correct mean, but the variances are off. I have tried it with other proposal distributions with roughly the same effect.

Thank you very much!

my code (sorry it's not properly vectorized):


clc
clear all
close all
numDim=2;
mu=[10 20];
sigma=[3 2.7;2.7 3];
sigma2=[.6 0;0 0.4];
numSamples=2000;
xnorm=zeros(numSamples,numDim);
mhsamplemult=10;
mhsamplebi=100;
pdf=@(x)((1/((2*pi)^(numDim/2)*det(sigma)^(1/2)))*exp(-(1/2)*(x'-mu')'*((sigma)\(x'-mu'))));
proppdf=[];
proprnd=@(x)(randn(size(x))*sigma2+x);
sym=1;

xvalstemp=mhsample(ones(1,numDim),numSamples,'pdf',pdf,'proprnd',proprnd,'proppdf',proppdf,'symmetric',sym,'thin',mhsamplemult,'burnin',mhsamplebi);
xmh=xvalstemp;

for ii=1:numSamples
xnew=mu+randn(1,numDim)*sigma;
xnorm(ii,:)=xnew;
end

for ii=1:numDim
figure
subplot(2,1,1)
hist(xnorm(:,ii),30)
title(['Distribution generated by randn in the ' num2str(ii) 'th variable'])
subplot(2,1,2)
hist(xmh(:,ii),30)
title(['Distribution generated by mhsample in the ' num2str(ii) 'th variable'])
end
if numDim==2
figure
subplot(2,1,1)
plot(xnorm(:,1),xnorm(:,2),'x')
title(['Distribution generated by randn'])
xlabel('x')
ylabel('y')
subplot(2,1,2)
plot(xmh(:,1),xmh(:,2),'x')
title(['Distribution generated by mhsample'])
xlabel('x')
ylabel('y')
end