% tutorial_mle.m
% Code from Myung (2003) Maximum Likelihood Tutorial
% Note: 'fmincon' function requires the MATLAB Optimization Toolbox
%
% Output:
%     POWER     EXPO
% w1  0.95312   1.0701  (0.88755 = r-squared POWER)
% w2  0.49793   0.13082 (0.96348 = r-squared EXPO)
% lik  313.37    305.31  exit1=5         exit2=4
%
% This is the main program that finds MLE estimates. Given a model, it
% takes sample size (n), time intervals (t) and observed proportion correct
% (y) as inputs. It returns the parameter values that maximize the log-
% likelihood function
clear all;
rand('seed',0);
global n t x; % define global variables
opts = optimset ('DerivativeCheck','off','Display','notify','TolX',1e-6,'TolFun',1e-6,...
   'Diagnostics','off','MaxIter',200,'LargeScale','off','Algorithm','active-set');
% option settings for optimization algorithm
n = 100 ;% number of independent Bernoulli trials (i.e., sample size)
t = [1 3 6 9 12 18]'; % time intervals as a column vector
y = [.94 .77 .40 .26 .24 .16]'; % observed proportion correct as a column vector
x = n*y; % number of correct responses

init_w = rand(2,1);   % starting parameter values
low_w = zeros(2,1);   % parameter lower bounds
up_w = 100*ones(2,1); % parameter upper bounds

display('Searching for maximimum likelihood values.  Please wait....');

while 1,
   [w1,lik1,exit1] = fmincon ('power_mle',init_w,[],[],[],[],low_w,up_w,[],opts);
   % optimization for power model that minimizes minus log-likelihood
   %(note that minimization of minus log-likelihood is equivalent to maximization 
   % of log-likelihood)
   %
   % w1: MLE parameter estimates
   % lik1: maximized log-likelihood value
   % exit1: optimization has converged if exit1 > 0 or not otherwise
   [w2,lik2,exit2] = fmincon('expo_mle',init_w,[],[],[],[],low_w,up_w,[],opts);
   % optimization for exponential model that minimizes minus log-likelihood
   prd1 = w1(1,1)*t.^(-w1(2,1)); % best fit prediction by power model
   r2(1,1) = 1-sum((prd1-y).^2)/sum((y-mean(y)).^2); % r^2 for power model
   prd2 = w2(1,1)*exp(-w2(2,1)*t); % best fit prediction by exponential model
   r2(2,1) = 1-sum((prd2-y).^2)/sum((y-mean(y)).^2); %r^2 for exponential model

   if sum(r2>0) == 2
      break; %exit infinite loop when both r-squared values exceed zero
   else
      init_w = rand(2,1); %otherwise reseed starting values and try again
   end;
end; %end(while)

format long;
display(' ');
display('Results:');
disp(num2str([w1 w2 r2],5));% display results
disp(num2str([lik1 lik2 exit1 exit2],5)); % display results
