[Octave-patch-tracker] [patch #10016] [octave forge] (statistics) mhsample implementation (over from bug #59924)

[anonymous]

Follow-up Comment #6, patch #10016 (project octave):

I decided to finish the m file implementation with the nchains propertie. I
made guesses on how I think Matlab would take inputs and present the output.
Unfortunately, Matlab's documentation does not give any examples using more
than one chain. The original submission had no issues with all of Matlab's
examples.

Metropolis-Hastings sampling is a method to generate random numbers from a
probability distribution. It has the following advantages over other sampling
methods:
* The distribution does not need to be normalized.
* It is particularly efficient in high dimensions where most of the pdf is
near zero (no curse of dimensionality), selection rejection sampling will
result in most values being rejected in this case
Its disadvantages are:
* Points next to each other in a chain are correlated.
* If the initial point is not in a location of high probability the initial
samples are not useful.

I have included below an example of sampling from a distribution made from
polynomials. I compare the sample mean, standard deviation along the axis
direction (diagonal), and the covariance (off-diagonal) with what I calculated
algebraically. This example demonstrated all of the inputs to the method and I
have commented on what I am unsure about. In the file, there are also two
demos. The first one samples a normal distribution and uses the sampled points
to integrate the top half of a sphere. The second one samples a single
variable normal distribution to find the normalization constant of a truncated
normal distribution. I have tested all of Matlab's examples.


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);
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]);
avg=mean(mxy,3)
fprintf("Diagonal of covariance answer: [%f %f]
\n",[(1/3478629)*sqrt(1110349)*sqrt(3478629)
(1/1159543)*sqrt(384653)*sqrt(1159543)]);
diagstd=mean(mean(sm.^2,1),3).^.5
fprintf("Off diagonal of covariance answer: %f \n",[-11346/165649]);
covar=mean(mean(prod(sm,2),1),3)
disp("Normally used to tune the method depends on proppdf and proprnd as well
as pdf");
avgacc=mean(accept)%Normally used for tuning the method


If the implementation lines up with Matlab I will provide the help for the
file and add the spaces to make it compliment with the Octave coding
guidelines. As the c++ version did not significantly speed up the function (a
few seconds faster for a 20-second runtime in pure Octave) I'm not sure if I
should update the c++ version.

The new version is attached.

(file #50790)
    _______________________________________________________

Additional Item Attachment:

File name: mhsample.m                     Size:4 KB
    <https://file.savannah.gnu.org/file/mhsample.m?file_id=50790>



    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/patch/?10016>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/