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