fixed typos and completed skeleton

Former-commit-id: 26ba426

Author: compops

matlab-seminar/complete/example1_iid.m

@@ -48,7 +48,7 @@
 
 hold on; 
 grid = 5:1:45;
-plot( grid, T .* exp( poisson( grid, thhat ) ), 'r', 'LineWidth', 2 ); 
+plot( grid, T .* exp( dpoisson( grid, thhat ) ), 'r', 'LineWidth', 2 ); 
 hold off;
 
 % Plot the parameter posterior estimate

matlab-seminar/complete/mh.m

@@ -19,16 +19,16 @@
 
 function[ theta ] = mh( observations, initialTheta, nIterations, Sigma )
 
- % Initalise variables
+ % Initialise variables
   theta         = zeros( nIterations, length(initialTheta) ); 
-  loglikelihood = zeros( nIterations, 1 );
+  logposterior  = zeros( nIterations, 1 );
   accept        = zeros( nIterations, 1 );
 
   % Set the initial parameter
   theta(1,:) = initialTheta;  
 
-  % Compute the initial log-likelihood
-  loglikelihood(1) = sum( dpoisson(observations, theta(1,:)) );
+  % Compute the initial log-posterior
+  logposterior(1) = sum( dpoisson(observations, theta(1,:)) );
 
   %=====================================================================
   % Run main loop
@@ -38,15 +38,15 @@
     % Propose a new parameter
     theta_proposed = mvnrnd( theta(kk-1,:), Sigma );
 
-    % Estimate the log-likelihood if the intensity is positive
+    % Compute the log-posterior if the intensity is positive
     if ( theta_proposed(1) > 0.0 )
-        loglikelihood_proposed = sum( dpoisson(observations, theta_proposed) );
+        logposterior_proposed = sum( dpoisson(observations, theta_proposed) );
     else
-        loglikelihood_proposed = -inf;
+        logposterior_proposed = -inf;
     end
 
     % Compute the acceptance probability
-    aprob = exp( loglikelihood_proposed - loglikelihood(kk-1) );
+    aprob = exp( logposterior_proposed - logposterior(kk-1) );
 
     % Generate uniform random variable in U[0,1]
     u = unifrnd(0,1);
@@ -56,13 +56,13 @@
 
       % Accept the parameter
       theta(kk,:)        = theta_proposed;
-      loglikelihood(kk)  = loglikelihood_proposed;
+      logposterior(kk)   = logposterior_proposed;
       accept(kk)         = 1.0;
 
     else
       % Reject the parameter
       theta(kk,:)        = theta(kk-1,:);
-      loglikelihood(kk)  = loglikelihood(kk-1);
+      logposterior(kk)   = logposterior(kk-1);
       accept(kk)         = 0.0;
     end
 

matlab-seminar/complete/pmh.m

@@ -21,16 +21,16 @@
 
 function[ theta ] = pmh( observations, initialTheta, nParticles, nIterations, Sigma )
 
- % Initalise variables
+ % Initialise variables
   theta          = zeros( nIterations, length(initialTheta) );
-  loglikelihood  = zeros( nIterations, 1 );
+  logposterior   = zeros( nIterations, 1 );
   accept         = zeros( nIterations, 1 );
 
   % Set the initial parameter
   theta(1,:) = initialTheta;  
 
-  % Run the initial particle filter to estimate the log-likelihood
-  [ ~ , loglikelihood(1) ] = pf( observations, theta(1,:), nParticles );
+  % Run the initial particle filter to estimate the log-posterior
+  [ ~ , logposterior(1) ] = pf( observations, theta(1,:), nParticles );
 
   %=====================================================================
   % Run main loop
@@ -40,15 +40,15 @@
     % Propose a new parameter
     theta_proposed = mvnrnd( theta(kk-1,:), Sigma );
 
-    % Estimate the log-likelihood (don't run if unstable system)
-    if ( abs( theta_proposed(1) ) < 1.0 )
-      [ ~, loglikelihood_proposed ] = pf( observations, theta_proposed, nParticles );
+    % Estimate the log-posterior (don't run if unstable system)
+    if ( ( abs( theta_proposed(1) ) < 1.0 ) && ( theta_proposed(2) > 0 ) && ( theta_proposed(3) > 0) )
+      [ ~, logposterior_proposed ] = pf( observations, theta_proposed, nParticles );
     else
-        loglikelihood_proposed = -inf;
+        logposterior_proposed = -inf;
     end
 
     % Compute the acceptance probability
-    aprob = exp( loglikelihood_proposed - loglikelihood(kk-1) );
+    aprob = exp( logposterior_proposed - logposterior(kk-1) );
 
     % Generate uniform random variable in U[0,1]
     u = unifrnd(0,1);
@@ -58,13 +58,13 @@
 
       % Accept the parameter
       theta(kk,:)       = theta_proposed;
-      loglikelihood(kk) = loglikelihood_proposed;
+      logposterior(kk)  = logposterior_proposed;
       accept(kk) = 1.0;
 
     else
       % Reject the parameter
       theta(kk,:)       = theta(kk-1,:);
-      loglikelihood(kk) = loglikelihood(kk-1);
+      logposterior(kk)  = logposterior(kk-1);
       accept(kk)        = 0.0;
     end
 

matlab-seminar/skeleton/dpoisson.m

@@ -0,0 +1,23 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Computes the log-PDF of Poisson distribution
+%
+% Inputs:
+% x:   Obserations
+% par: Parameters of the PDF (location, scale and shape)
+%
+% Outputs:
+% out: Logarithm of the PDF evaluated at each x.
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function[out] = poisson(x, par)
+    part1 = -log( factorial( x ) );
+    part2 = x * log( par );
+    part3 = -par;
+    out = part1 + part2 + part3;
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% End of file
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

matlab-seminar/skeleton/dt.m

@@ -17,18 +17,15 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 % The inital guess of the parameter
-initPar  = [ 0.5 0.5 15 ];
+initPar  = mean(y);
 
 % The length of the burn-in and the no. iterations of MH algorithm 
 % ( nBurnIn < nIterations )
 nBurnIn     = 10000;
 nIterations = 50000;
 
 % The covariance matrix in the random walk proposal
-Sigma = [ [ 121.3141    1.4331    3.1442];
-          [   1.4331    0.4692    0.0443];
-          [   3.1442    0.0443    0.3539] ];
-Sigma = 0.8 * Sigma;
+Sigma = 0.17;
 
 % Run the MH algorithm
 th = mh( y, initPar, nIterations, Sigma );
@@ -44,42 +41,26 @@
 figure(2);
 
 % Plot the data and the fitted model
-subplot(4,2,[1 2]);
+subplot(2,2,[1 2]);
 hist( y, 25 )
 xlabel( 'time' ); 
 ylabel( 'no. earthquakes' );
 
 hold on; 
-grid = 5:0.1:45;
-plot( grid, T .* exp( dt( grid, thhat ) ), 'r', 'LineWidth', 2 ); 
+grid = 5:1:45;
+plot( grid, T .* exp( poisson( grid, thhat ) ), 'r', 'LineWidth', 2 ); 
 hold off;
 
 % Plot the parameter posterior estimate
 % Plot the trace of the Markov chain after burn-in
-subplot(4,2,3);
+subplot(2,2,3);
 hist( th( nBurnIn:nIterations, 1 ), floor( sqrt( nIterations - nBurnIn ) ) );
-xlabel('nu'); 
+xlabel('theta'); 
 ylabel('posterior density estimate');
 
-subplot(4,2,4);
+subplot(2,2,4);
 plot( nBurnIn:nIterations, th( nBurnIn:nIterations, 1 ) );
-xlabel('iteration'); ylabel('trace of nu');
-
-subplot(4,2,5);
-hist( th( nBurnIn:nIterations, 2 ), floor( sqrt( nIterations - nBurnIn ) ) );
-xlabel('mu'); ylabel('posterior density estimate');
-
-subplot(4,2,6);
-plot( nBurnIn:nIterations, th( nBurnIn:nIterations, 2 ) );
-xlabel('iteration'); ylabel('trace of mu');
-
-subplot(4,2,7);
-hist( th( nBurnIn:nIterations, 3 ), floor( sqrt( nIterations - nBurnIn ) ) );
-xlabel('sigma'); ylabel('posterior density estimate');
-
-subplot(4,2,8);
-plot( nBurnIn:nIterations, th( nBurnIn:nIterations, 3 ) );
-xlabel('iteration'); ylabel('trace of sigma');
+xlabel('iteration'); ylabel('trace of theta');
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % End of file

matlab-seminar/skeleton/example1_iid.m

@@ -19,34 +19,54 @@
 
 function[ theta ] = mh( observations, initialTheta, nIterations, Sigma )
 
- % Initalise variables
+ % Initialise variables
   theta         = zeros( nIterations, length(initialTheta) ); 
-  loglikelihood = zeros( nIterations, 1 );
+  logposterior  = zeros( nIterations, 1 );
   accept        = zeros( nIterations, 1 );
 
   % Set the initial parameter
   theta(1,:) = initialTheta;  
 
-  % Compute the initial log-likelihood
-  loglikelihood(1) = sum( dt(observations, theta(1,:)) );
+  % Compute the initial log-posterior
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  %%%%%%%%%%%%%%%%% ADD CODE HERE %%%%%%%%%%%%%%%%%%%%%
+  % Hint: dpossion and observations                   %
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  logposterior(1) = ...
 
   %=====================================================================
   % Run main loop
   %=====================================================================
   for kk = 2:nIterations
 
     % Propose a new parameter
-    theta_proposed = mvnrnd( theta(kk-1,:), Sigma );
 
-    % Estimate the log-likelihood if DOF and sigma is positive
-    if ( ( theta_proposed(1) > 0.0 ) && ( theta_proposed(3) > 0.0 ) )
-        loglikelihood_proposed = sum( dt(observations, theta_proposed) );
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    %%%%%%%%%%%%%%%%% ADD CODE HERE %%%%%%%%%%%%%%%%%%%%%
+    % Hint: mvnrnd and Sigma                            %
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    theta_proposed = ...
+    
+    % Compute the log-posterior if the intensity is positive
+    if ( theta_proposed(1) > 0.0 )
+        
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        %%%%%%%%%%%%%%%%% ADD CODE HERE %%%%%%%%%%%%%%%%%%%%%
+        % Hint: dpossion and observations                   %
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        logposterior_proposed = ...
+          
     else
-        loglikelihood_proposed = -inf;
+        logposterior_proposed = -inf;
     end
 
     % Compute the acceptance probability
-    aprob = exp( loglikelihood_proposed - loglikelihood(kk-1) );
+    
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    %%%%%%%%%%%%%%%%% ADD CODE HERE %%%%%%%%%%%%%%%%%%%%%
+    % Hint: exp, logposterior_proposed and logposterior
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    aprob = ...
 
     % Generate uniform random variable in U[0,1]
     u = unifrnd(0,1);
@@ -56,13 +76,13 @@
 
       % Accept the parameter
       theta(kk,:)        = theta_proposed;
-      loglikelihood(kk)  = loglikelihood_proposed;
+      logposterior(kk)   = logposterior_proposed;
       accept(kk)         = 1.0;
 
     else
       % Reject the parameter
       theta(kk,:)        = theta(kk-1,:);
-      loglikelihood(kk)  = loglikelihood(kk-1);
+      logposterior(kk)   = logposterior(kk-1);
       accept(kk)         = 0.0;
     end