function fit = fit_meta_d_MLE(nR_S1, nR_S2, s, fncdf, fninv)

% fit = fit_meta_d_MLE(nR_S1, nR_S2, s, fncdf, fninv)
%
% Given data from an experiment where an observer discriminates between two
% stimulus alternatives on every trial and provides confidence ratings,
% provides a type 2 SDT analysis of the data.
%
% INPUTS
%
% * nR_S1, nR_S2
% these are vectors containing the total number of responses in
% each response category, conditional on presentation of S1 and S2.
%
% e.g. if nR_S1 = [100 50 20 10 5 1], then when stimulus S1 was
% presented, the subject had the following response counts:
% responded S1, rating=3 : 100 times
% responded S1, rating=2 : 50 times
% responded S1, rating=1 : 20 times
% responded S2, rating=1 : 10 times
% responded S2, rating=2 : 5 times
% responded S2, rating=3 : 1 time
%
% The ordering of response / rating counts for S2 should be the same as it
% is for S1. e.g. if nR_S2 = [3 7 8 12 27 89], then when stimulus S2 was
% presented, the subject had the following response counts:
% responded S1, rating=3 : 3 times
% responded S1, rating=2 : 7 times
% responded S1, rating=1 : 8 times
% responded S2, rating=1 : 12 times
% responded S2, rating=2 : 27 times
% responded S2, rating=3 : 89 times
%
% * s
% this is the ratio of standard deviations for type 1 distributions, i.e.
%
% s = sd(S1) / sd(S2)
%
% if not specified, s is set to a default value of 1.
% the function SDT_MLE_fit available at
% http://www.columbia.edu/~bsm2105/type2sdt/
% can be used to get an MLE estimate of s using Matlab's optimization
% toolbox
%
% * fncdf
% a function handle for the CDF of the type 1 distribution.
% if not specified, fncdf defaults to @normcdf (i.e. CDF for normal
% distribution)
%
% * fninv
% a function handle for the inverse CDF of the type 1 distribution.
% if not specified, fninv defaults to @norminv
%
% OUTPUT
%
% Output is packaged in the struct "fit." 
% In the following, let S1 and S2 represent the distributions of evidence 
% generated by stimulus classes S1 and S2.
% Then the fields of "fit" are as follows:
% 
% fit.da        = mean(S2) - mean(S1), in room-mean-square(sd(S1),sd(S2)) units
% fit.s         = sd(S1) / sd(S2)
% fit.meta_da   = meta-d' in RMS units
% fit.M_diff    = meta_da - da
% fit.M_ratio   = meta_da / da
% fit.meta_ca   = type 1 criterion for meta-d' fit, RMS units
% fit.t2ca_rS1  = type 2 criteria of "S1" responses for meta-d' fit, RMS units
% fit.t2ca_rS2  = type 2 criteria of "S2" responses for meta-d' fit, RMS units
%
% fit.S1units   = contains same parameters in sd(S1) units.
%                 these may be of use since the data-fitting is conducted  
%                 using parameters specified in sd(S1) units.
% 
% fit.logL          = log likelihood of the data fit
%
% fit.est_HR2_rS1  = estimated (from meta-d' fit) type 2 hit rates for S1 responses
% fit.obs_HR2_rS1  = actual type 2 hit rates for S1 responses
% fit.est_FAR2_rS1 = estimated type 2 false alarm rates for S1 responses
% fit.obs_FAR2_rS1 = actual type 2 false alarm rates for S1 responses
% 
% fit.est_HR2_rS2  = estimated type 2 hit rates for S2 responses
% fit.obs_HR2_rS2  = actual type 2 hit rates for S2 responses
% fit.est_FAR2_rS2 = estimated type 2 false alarm rates for S2 responses
% fit.obs_FAR2_rS2 = actual type 2 false alarm rates for S2 responses
%
% If there are N ratings, then there will be N-1 type 2 hit rates and false
% alarm rates. 

% 9/7/10 - bm - wrote it


%% parse inputs

if ~mod(length(nR_S1),2)==0, error('input arrays must have an even number of elements'); end
if length(nR_S1)~=length(nR_S2), error('input arrays must have the same number of elements'); end

if ~exist('s','var') || isempty(s)
    s = 1;
end

if ~exist('fncdf','var') || isempty(fncdf)
    fncdf = @normcdf;
end

if ~exist('fninv','var') || isempty(fninv)
    fninv = @norminv;
end

nRatings = length(nR_S1) / 2;
nCriteria = 2*nRatings - 1;


%% set up constraints for MLE estimation

% parameters
% meta-d' - 1
% t2c     - nCriteria-1

A = [];
b = [];

% constrain type 2 criteria values,
% such that t2c(i) is always <= t2c(i+1)
% want t2c(i)   <= t2c(i+1) 
% -->  t2c(i+1) >= c(i) + 1e-5 (i.e. very small deviation from equality) 
% -->  t2c(i) - t2c(i+1) <= -1e-5 
for i = 2 : nCriteria-1

    A(end+1,[i i+1]) = [1 -1];
    b(end+1) = -1e-5;

end

% lower bounds on parameters
LB = [];
LB = [LB -10];                         % meta-d'
LB = [LB -20*ones(1,(nCriteria-1)/2)]; % criteria lower than t1c
LB = [LB zeros(1,(nCriteria-1)/2)];    % criteria higher than t1c

% upper bounds on parameters
UB = [];
UB = [UB 10];                          % meta-d'
UB = [UB zeros(1,(nCriteria-1)/2)];    % criteria lower than t1c
UB = [UB 20*ones(1,(nCriteria-1)/2)];  % criteria higher than t1c


%% select constant criterion type

constant_criterion = 'meta_d1 * (t1c1 / d1)'; % relative criterion


%% set up initial guess at parameter values

ratingHR  = [];
ratingFAR = [];
for c = 2:nRatings*2
    ratingHR(end+1) = sum(nR_S2(c:end)) / sum(nR_S2);
    ratingFAR(end+1) = sum(nR_S1(c:end)) / sum(nR_S1);
end

t1_index = nRatings;
t2_index = setdiff(1:2*nRatings-1, t1_index);

d1 = (1/s) * fninv( ratingHR(t1_index) ) - fninv( ratingFAR(t1_index) );
meta_d1 = d1;

c1 = (-1/(1+s)) * ( fninv( ratingHR ) + fninv( ratingFAR ) );
t1c1 = c1(t1_index);
t2c1 = c1(t2_index);

guess = [meta_d1 t2c1 - eval(constant_criterion)];



%% find the best fit for type 2 hits and FAs

% save fit_meta_d_MLE.mat nR_S1 nR_S2 t2FAR_rS2 t2HR_rS2 t2FAR_rS1 t2HR_rS1 nRatings t1c1 s d1 fncdf constant_criterion
save fit_meta_d_MLE.mat nR_S1 nR_S2 nRatings d1 t1c1 s constant_criterion fncdf fninv

op = optimset(@fmincon);
op = optimset(op,'MaxFunEvals',100000);

[x f] = fmincon(@fit_meta_d_logL,guess,A,b,[],[],LB,UB,[],op);

meta_d1  = x(1);
t2c1     = x(2:end) + eval(constant_criterion);
logL     = -f;


%% data is fit, now to package it...

%% find observed t2FAR and t2HR 

% I_nR and C_nR are rating trial counts for incorrect and correct trials
% element i corresponds to # (in)correct w/ rating i
I_nR_rS2 = nR_S1(nRatings+1:end);
I_nR_rS1 = nR_S2(nRatings:-1:1);

C_nR_rS2 = nR_S2(nRatings+1:end);
C_nR_rS1 = nR_S1(nRatings:-1:1);

for i = 2:nRatings
    obs_FAR2_rS2(i-1) = sum( I_nR_rS2(i:end) ) / sum(I_nR_rS2);
    obs_HR2_rS2(i-1)  = sum( C_nR_rS2(i:end) ) / sum(C_nR_rS2);
    
    obs_FAR2_rS1(i-1) = sum( I_nR_rS1(i:end) ) / sum(I_nR_rS1);
    obs_HR2_rS1(i-1)  = sum( C_nR_rS1(i:end) ) / sum(C_nR_rS1);    
end


%% find estimated t2FAR and t2HR

S1mu = -meta_d1/2; S1sd = 1;
S2mu =  meta_d1/2; S2sd = S1sd/s;

mt1c1 = eval(constant_criterion);

C_area_rS2 = 1-fncdf(mt1c1,S2mu,S2sd);
I_area_rS2 = 1-fncdf(mt1c1,S1mu,S1sd);

C_area_rS1 = fncdf(mt1c1,S1mu,S1sd);
I_area_rS1 = fncdf(mt1c1,S2mu,S2sd);

for i=1:nRatings-1
    
    t2c1_lower = t2c1(nRatings-i);
    t2c1_upper = t2c1(nRatings-1+i);
        
    I_FAR_area_rS2 = 1-fncdf(t2c1_upper,S1mu,S1sd);
    C_HR_area_rS2  = 1-fncdf(t2c1_upper,S2mu,S2sd);

    I_FAR_area_rS1 = fncdf(t2c1_lower,S2mu,S2sd);
    C_HR_area_rS1  = fncdf(t2c1_lower,S1mu,S1sd);
    
    
    est_FAR2_rS2(i) = I_FAR_area_rS2 / I_area_rS2;
    est_HR2_rS2(i)  = C_HR_area_rS2 / C_area_rS2;

    est_FAR2_rS1(i) = I_FAR_area_rS1 / I_area_rS1;
    est_HR2_rS1(i)  = C_HR_area_rS1 / C_area_rS1;
    
end


%% package output

fit.da        = sqrt(2/(1+s^2)) * s * d1;
fit.s         = s;
fit.meta_da   = sqrt(2/(1+s^2)) * s * meta_d1;
fit.M_diff    = fit.meta_da - fit.da;
fit.M_ratio   = fit.meta_da / fit.da;
fit.meta_ca   = ( sqrt(2).*s ./ sqrt(1+s.^2) ) .* t1c1;
t2ca          = ( sqrt(2).*s ./ sqrt(1+s.^2) ) .* t2c1;
fit.t2ca_rS1  = t2ca(1:nRatings-1);
fit.t2ca_rS2  = t2ca(nRatings:end);
fit.S1units.d1        = d1;
fit.S1units.meta_d1   = meta_d1;
fit.S1units.s         = s;
fit.S1units.meta_c1   = t1c1;
fit.S1units.t2c1_rS1  = t2c1(1:nRatings-1);
fit.S1units.t2c1_rS2  = t2c1(nRatings:end);

fit.logL           = logL;

fit.est_HR2_rS1  = est_HR2_rS1;
fit.obs_HR2_rS1  = obs_HR2_rS1;

fit.est_FAR2_rS1 = est_FAR2_rS1;
fit.obs_FAR2_rS1 = obs_FAR2_rS1;

fit.est_HR2_rS2  = est_HR2_rS2;
fit.obs_HR2_rS2  = obs_HR2_rS2;

fit.est_FAR2_rS2 = est_FAR2_rS2;
fit.obs_FAR2_rS2 = obs_FAR2_rS2;


%% clean up
delete fit_meta_d_MLE.mat

end


%% function to find the likelihood of parameter values, given observed data
function logL = fit_meta_d_logL(parameters)

% set up parameters
meta_d1  = parameters(1);
t2c1     = parameters(2:end);

% loads:
% nR_S1 nR_S2 nRatings d1 t1c1 s constant_criterion fncdf fninv
load fit_meta_d_MLE.mat


% define mean and SD of S1 and S2 distributions
S1mu = -meta_d1/2; S1sd = 1;
S2mu =  meta_d1/2; S2sd = S1sd/s;


% adjust so that the type 1 criterion is set at 0
% (this is just to work with optimization toolbox constraints...
%  to simplify defining the upper and lower bounds of type 2 criteria)
S1mu = S1mu - eval(constant_criterion);
S2mu = S2mu - eval(constant_criterion);

t1c1 = 0;



%%% set up MLE analysis

% get type 2 response counts
for i = 1:nRatings
    
    % S1 responses
    nC_rS1(i) = nR_S1(i);
    nI_rS1(i) = nR_S2(i);
    
    % S2 responses
    nC_rS2(i) = nR_S2(nRatings+i);
    nI_rS2(i) = nR_S1(nRatings+i);
    
end

% get type 2 probabilities
C_area_rS1 = fncdf(t1c1,S1mu,S1sd);
I_area_rS1 = fncdf(t1c1,S2mu,S2sd);

C_area_rS2 = 1-fncdf(t1c1,S2mu,S2sd);
I_area_rS2 = 1-fncdf(t1c1,S1mu,S1sd);

t2c1x = [-Inf t2c1(1:nRatings-1) t1c1 t2c1(nRatings:end) Inf];

for i = 1:nRatings
    prC_rS1(i) = ( fncdf(t2c1x(i+1),S1mu,S1sd) - fncdf(t2c1x(i),S1mu,S1sd) ) / C_area_rS1;
    prI_rS1(i) = ( fncdf(t2c1x(i+1),S2mu,S2sd) - fncdf(t2c1x(i),S2mu,S2sd) ) / I_area_rS1;
    
    prC_rS2(i) = ( (1-fncdf(t2c1x(nRatings+i),S2mu,S2sd)) - (1-fncdf(t2c1x(nRatings+i+1),S2mu,S2sd)) ) / C_area_rS2;
    prI_rS2(i) = ( (1-fncdf(t2c1x(nRatings+i),S1mu,S1sd)) - (1-fncdf(t2c1x(nRatings+i+1),S1mu,S1sd)) ) / I_area_rS2;
end
    

% calculate logL
logL = 0;
for i = 1:nRatings
    
    logL = logL + nC_rS1(i)*log(prC_rS1(i)) + nI_rS1(i)*log(prI_rS1(i)) + ...
                  nC_rS2(i)*log(prC_rS2(i)) + nI_rS2(i)*log(prI_rS2(i));
              
end

if isnan(logL), logL=-Inf; end

logL = -logL;
  
end