nchain = 1000; % Number of chains
% Not sure how the starting point is input all three of the below work for my implementation
% I allow different starting points, I hope Matlab does as well
%start=rand(nchain,2);
start = rand (1, 2, nchain);
%start=rand(1,2);%all chains start from same point
nsamples = 10000; % Length of chain
pdf = @(x)(-(x(:, 1) - 1 / 2) .^ 2 .* (x(:, 1) - 1) .* (x(:, 1) + 1) - ...
x(:, 1) .* x(:, 2) - (x(:, 2) + 1/3) .^ 2 .* (x(:, 2) + 1) .* ...
(x(:, 2) - 1) + 1) .* (x(:, 1) >= -1 & x(:, 1) <= 1 & ...
x(:, 2) >= -1 & x(:, 2) <= 1); % * 135 / 814; % The distribution we wish to sample from, does not need to be normalized. Also assuming vectorized function evaluation have different points down and different variables across
delta = .5;
proppdf = @(x, y) prod (unifpdf (y - x, -delta, delta), 2); % pdf to choose next point in chain
proprnd = @(x) x + rand (size (x)) * 2 * delta - delta; %function to draw random number from above pdf
sym = true; % if proppdf is a symmetric distribution proppdf not used
K = 0; % number of samples to discard at the beginning
m = 2; %disgard (m-1) every m samples
% not sure of the output sizes of smpl and accept. I have values of a chain
% dimension 1, variables dimension 2 and different chains dimension 3
[smpl, accept] = mhsample (start, nsamples, 'pdf', pdf, 'proppdf', proppdf, ...
'proprnd',proprnd, 'symmetric', sym, 'burnin', K, ...
'thin', m, 'nchain', nchain);
{The logical indices in position 1 contain a true value outside of the array bounds.
Error in mhsample (line 178)
x0(acc,:) = y(acc,:); % preserves x's shape.
}
[smpl,accept]=mhsample(start,nsamples,'pdf',pdf,'proppdf',proppdf,'proprnd',proprnd,'symmetric',sym,'burnin',K,'thin',m,'nchain',nchain);
{The logical indices in position 1 contain a true value outside of the array bounds.
Error in mhsample (line 178)
x0(acc,:) = y(acc,:); % preserves x's shape.
}
nchain = 1000; % Number of chains
% Not sure how the starting point is input all three of the below work for my implementation
% I allow different starting points, I hope Matlab does as well
start=rand(nchain,2);
%start = rand (1, 2, nchain);
%start=rand(1,2);%all chains start from same point
nsamples = 10000; % Length of chain
% The distribution we wish to sample from, does not need to be normalized.
% Also assuming vectorized function evaluation have different points down
% and different variables across
pdf = @(x) (-(x(:, 1) - 1 / 2) .^ 2 .* (x(:, 1) - 1) .* (x(:, 1) + 1) - ...
x(:, 1) .* x(:, 2) - (x(:, 2) + 1/3) .^ 2 .* (x(:, 2) + 1) .* ...
(x(:, 2) - 1) + 1) .* (x(:, 1) >= -1 & x(:, 1) <= 1 & ...
x(:, 2) >= -1 & x(:, 2) <= 1); % * 135 / 814;
delta = .5;
proppdf = @(x, y) prod (unifpdf (y - x, -delta, delta), 2); % pdf to choose next point in chain
proprnd = @(x) x + rand (size (x)) * 2 * delta - delta; %function to draw random number from above pdf
sym = true; % if proppdf is a symmetric distribution proppdf not used
K = 0; % number of samples to discard at the beginning
m = 2; %disgard (m-1) every m samples
% not sure of the output sizes of smpl and accept. I have values of a chain
% dimension 1, variables dimension 2 and different chains dimension 3
[smpl, accept] = mhsample (start, nsamples, 'pdf', pdf, 'proppdf', proppdf, ...
'proprnd',proprnd, 'symmetric', sym, 'burnin', K, ...
'thin', m, 'nchain', nchain);
close all;
hold on;
[xx yy] = meshgrid (linspace (-1, 1, 25));
mesh (xx, yy, 135 / 814 * reshape (pdf ([xx(:), yy(:)]), size (xx)), 'facecolor', 'none');
plot (smpl(:, 1, 1), smpl(:, 2, 1), '.');
hold off;
mxy = mean (smpl, 1);
sm = (smpl - mxy);
fprintf ("Mean answer: [%f %f] \n", [-36/407 24/407]);
Mean answer: [-0.088452 0.058968]
avg = mean (mxy, 3)
avg =
-0.0885 0.0588
fprintf ("Diagonal of covariance answer: [%f %f] \n", ...
[((1 / 3478629) * sqrt (1110349) * sqrt (3478629)) ...
((1 / 1159543) * sqrt (384653) * sqrt (1159543))]);
Diagonal of covariance answer: [0.564970 0.575958]
diagstd = mean (mean (sm .^ 2, 1), 3) .^ .5
diagstd =
0.5645 0.5758
fprintf ("Off diagonal of covariance answer: %f \n", [-11346/165649]);
Off diagonal of covariance answer: -0.068494
covar = mean (mean (prod (sm, 2), 1), 3)
covar =
-0.0685
disp ("Normally used to tune the method depends on proppdf and proprnd as well as pdf");
Normally used to tune the method depends on proppdf and proprnd as well as pdf
avgacc = mean (accept)
avgacc =
0.7083
diary off