
clear
clf

% Monte Carlo sampling of a probability density function
% using the simulated annealing approach

% Heiner Igel, 2002

% 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, gravity inversion, etc. )

% Some global parameters
np=1000000;				% number of tries to take
ns=100;				% number of models to keep
T0=1.;
alpha=0.995;
T=T0;

% 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)));
ppdf=pdf*0+1;
ppdf=ppdf/(sum(sum(ppdf)));


T=input(' T ');

% calculate the a posteriori pdf as afunction of T

for i=1:n,
   for j=1:n,
      S(i,j)=-T0*log(pdf(i,j)/ppdf(i,j));
      smT(i,j)=ppdf(i,j)*exp( -S(i,j)/T);
   end
end


pdf=pdf/max(max(pdf));
smT=smT/max(max(smT));
S=S/max(max(S));

subplot(131),surf(pdf),axis([1 n 1 n 0 1]), shading interp, light, axis square
subplot(132),surf(smT),axis([1 n 1 n 0 1]), shading interp, light, axis square
subplot(133),surf(S),axis([1 n 1 n 0 1]), shading interp, light, axis square



