clear
clf

% Monte Carlo sampling of a probability density function

% for simplicity we use the "peaks" function. The corresponding
% section can simply be replaced with the probability from a
% proper inverse problem (e.g. hypocenter location)

% modified to look only in the near neighborhood

% Some global parameters
%np=1000000;				% number of tries to take
ns=100;				% number of samples to keep 
neigh=1;					% look only within 5% of the total size of the model space

% Let's define a pdf

n=100;				% number of dimension in x and y
pdf=peaks(n);
pdf=sqrt(pdf.*pdf);
pdf=pdf/(sum(sum(pdf)));

surf(pdf)

% Random walk

% find an initial vector x

xcur=round([rand rand]*n);

% Let s get moving

is=1;
np=0;
xa(1,1)=xcur(1);
xa(1,2)=xcur(2);

imagesc(pdf'), axis off 
hold on

while is<=ns,

	np=np+1;   
   % make a random choice for the next move
   xnew=xcur+round(([rand rand]-.5)*n*neigh);
   xnew(find(xnew<=0))=xnew(find(xnew<=0))*0+1;
   xnew(find(xnew>n))=xnew(find(xnew>n))*0+n;
   
   % compare probabilities
   Pcur=pdf(xcur(1),xcur(2));
   Pnew=pdf(xnew(1),xnew(2));
   
   if Pnew >= Pcur, 
      xcur=xnew;
      is=is+1;
      np;
      if (xcur(1)==77 & xcur(2)==50), 
          is
          np
          break;end;
      xa(is,1)=xcur(1);
      xa(is,2)=xcur(2);
      Pa(is)=Pnew;
      disp(sprintf(' Made the %i-th move to [%i,%i] (Pnew>Pold) ', is,xcur(1),xcur(2))) 
      
      % display move graphically
      
      plot([xa(is-1,1) xa(is,1)],[xa(is-1,2) xa(is,2)],'k-')
      plot([xa(is-1,1) xa(is,1)],[xa(is-1,2) xa(is,2)],'k+')
      drawnow
      
   end
   
   if Pnew < Pcur, 
      P=Pnew/Pcur;
      test=rand;
      if test<=P, 
      xcur=xnew;
      is=is+1;
      if (xcur(1)==77 & xcur(2)==50), 
          is
          np
          break;end;
      xa(is,1)=xcur(1);
      xa(is,2)=xcur(2);
       Pa(is)=Pnew;
      disp(sprintf(' Made the %i-th move to [%i,%i] (Pnew<Pold) ', is,xcur(1),xcur(2))) 
      
      % display move graphically
      
      plot([xa(is-1,1) xa(is,1)],[xa(is-1,2) xa(is,2)],'k-')
      plot([xa(is-1,1) xa(is,1)],[xa(is-1,2) xa(is,2)],'k+')
      drawnow
      
     end
   end  
end
np

hold off

% Final sampling

%imagesc(pdf'), hold on
%plot(xa(1:ns,1),xa(1:ns,2),'w+')