changed to Poisson IID data

Former-commit-id: 08805b9

Author: compops

matlab-seminar/complete/example1_iid.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/complete/mh.m

@@ -28,7 +28,7 @@
   theta(1,:) = initialTheta;  
 
   % Compute the initial log-likelihood
-  loglikelihood(1) = sum( dt(observations, theta(1,:)) );
+  loglikelihood(1) = sum( dpoisson(observations, theta(1,:)) );
 
   %=====================================================================
   % Run main loop
@@ -38,9 +38,9 @@
     % 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) );
+    % Estimate the log-likelihood if the intensity is positive
+    if ( theta_proposed(1) > 0.0 )
+        loglikelihood_proposed = sum( dpoisson(observations, theta_proposed) );
     else
         loglikelihood_proposed = -inf;
     end