%% Simple Metropolis algorithm example
%{
---------------------------------------------------------------------------
*Created by:                  Date:            Comment:
 Felipe Uribe-Castillo        July 2011        Final project algorithm  
*Hacked by:
Amy Graves                    March 2018       For Phys 114 Seminar
*Mail: 
 felipeuribecastillo@gmx.com
 abug1@swarthmore.edu
*
---------------------------------------------------------------------------
The Metropolis algorithm:
Obtains samples from an arbitrary probability distribution

Original references: 
-"The Monte Carlo method"
  Metropolis & Ulam (1949).
-"Equations of State Calculations by Fast Computing Machines"
  Metropolis, Rosenbluth, Rosenbluth, Teller & Teller (1953).
---------------------------------------------------------------------------
More recent references:
1."Markov chain Monte Carlo and Gibbs sampling"
   B.Walsh ----- Lecture notes for EEB 581, version april 2004.
   http://web.mit.edu/~wingated/www/introductions/mcmc-gibbs-intro.pdf
%}
clear; close all; home; format long g;

%% Initial data
% Probability Distribution information


beta = 2.0;                                             % A Target PDF parameter
pBoltz =  @(t) exp(-beta .* t);                              % A  Target "PDF" ... Boltzmann distribution

pUniform = @(t) t.^0 ;                                  % A  Target "PDF" ... uniform 
                                       
                                              
                                                        

nu           = 5;                                       % A Target PDF parameter
pGauss           = @(t) exp(- .5 * (t - nu).^2) ;                       % A  Target "PDF" ... Gaussian distribution



proposal_PDF = @(mu) unifrnd(0,10);                      % Proposal PDF
% Parameters
N        = 5000;              % Number of samples (iterations)
theta    = zeros(1,N);        % Samples drawn form the target "PDF"
theta(1) = 0.3;               % Initial state of the chain

%% Procedure

% Talk to user, and prompt for input
disp('Menu of PDFs:  Boltz with beta=2; uniform; Gaussian with mean=5, sigma = 1');
prompt = 'Choose your distribution 1:Boltz, 2: Uniform, 3: Gaussian \n' ;
d = input(prompt);
if(d == 1) p = pBoltz;
    elseif(d == 2) p = pUniform;
    elseif (d == 3) p = pGauss;
    else disp(' I did not understand  ... using Boltz') 
         p = pBoltz;
    end

    
    
for i = 1:N
   theta_ast = proposal_PDF([]); %currently, proposal_PDF is 10*rand
   % Sampling from the proposal PDF
   alpha     = p(theta_ast)/p(theta(i));  % Ratio of the density at theta_ast and theta points
   if rand <= min(alpha,1)
      % Accept the sample with prob = min(alpha,1)
      theta(i+1) = theta_ast;
   else
      % Reject the sample with prob = 1 - min(alpha,1)
      theta(i+1) = theta(i);
   end
end



% Samples and distributions
xx = eps:0.01:10;   % x-axis (Graphs)
figure;
% Histogram and target dist
subplot(2,1,2); 
[n1 x1] = hist(theta, ceil(sqrt(N)));
bar(x1,n1/(N*(x1(2)-x1(1))));   colormap summer;   hold on;   % Normalized histogram
plot(xx, p(xx)/trapz(xx,p(xx)), 'r-', 'LineWidth', 3);        % Normalized function
grid on;   xlim([0 max(theta)]);   
title('Distribution of samples', 'FontSize', 14);
ylabel('Probability density function', 'FontSize', 12);
% Samples
subplot(2,1,1);     
plot(theta, 1:N+1, 'b-');   xlim([0 max(theta)]);   ylim([0 N]);   grid on;
xlabel('Location', 'FontSize', 12);
ylabel('Iterations, N', 'FontSize', 12); 

%%End