function [analysis params] = SDT_MLE_fit(nR_S1, nR_S2)

% [analysis params] = SDT_MLE_fit(nR_S1, nR_S2)
% 
% Find a maximum-likelihood fit of the unequal variance signal detection 
% theory model to the input 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
%
% 
% OUTPUTS
% 
% * analysis
% This is a struct containing information about the data fitting procedure.
% Its fields are:
% 
% analysis.logL - the log likelihood of the fit
% analysis.k    - number of model parameters (d', s, and criteria)
% analysis.n    - total number of data points being estimated
% analysis.AIC  - Akaike Information Criterion, -2*logL + 2*k
% analysis.AICc - Small sample correction for AIC, -2*logL + (2*k*n)/(n-k-1)
% analysis.BIC  - Bayesian Information Criterion, -2*logL + k*log(n)
%
% * params
% This is a struct containing the MLE parameter values for the SDT model.
% In the following, let S1 and S2 represent the distributions of evidence 
% generated by stimulus classes S1 and S2.
% Then the params fields are:
%
% params.d1    - mean(S2) - mean(S1), in sd(S1) units
% params.c1    - decision criteria, in sd(S1) units
% params.s     - sd(S1) / sd(S2)
% params.da    - mean(S2) - mean(S1), in room-mean-square(sd(S1),sd(S2)) units
% params.ca    - decision criteria (in RMS units)
% params.guess - the initial guess for model parameter values 
%                (d1, s, and c1)


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


%% write the data for logL function

save current_SDTMLE_data.mat nR_S1 nR_S2 nCriteria nRatings


%% set up constraints
A = [];
b = [];

% constrain criteria values,
% such that c(i) is always <= c(i+1)
% want c(i)   <= c(i+1) 
% -->  c(i+1) >= c(i) + 1e-5 (i.e. very small deviation from equality) 
% -->  c(i) - c(i+1) <= -1e-5 
offset = 2;
for crit = 1 : nCriteria-1
    % c(crit) <= c(crit+1) --> c(crit) - c(crit+1) <= -1e-5
    A(end+1,[offset+crit offset+crit+1]) = [1 -1];
    b(end+1) = -1e-5;
end

% lower bounds on parameter values
LB = [];
LB = [LB -10];       % d1
LB = [LB 1/10];      % s
LB = [LB -20*ones(1,nCriteria)]; % c1

% upper bounds on parameter values
UB = [];
UB = [UB 10];        % d1
UB = [UB 10];        % s
UB = [UB 20*ones(1,nCriteria)];  % c1



%% 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

if isempty(ratingFAR) || all(ratingFAR == ratingFAR(1))
    s = 1;
else
    p = polyfit(norminv(ratingFAR), norminv(ratingHR), 1);
    s = p(1);
end

t1_index = nRatings;

d1 = (1/s) * norminv( ratingHR ) - norminv( ratingFAR );
% d1 = d1(t1_index);
d1 = mean(d1);

c1 = (-1/(1+s)) * ( norminv( ratingHR ) + norminv( ratingFAR ) );

guess = [d1 s c1];



%% do maximum likelihood estimation and package output

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

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

logL = -f;
k    = length(guess);
n    = length(nR_S1) + length(nR_S2);

analysis.logL = logL;
analysis.k    = k;
analysis.n    = n;
analysis.AIC  = -2*logL + 2*k;
analysis.AICc = -2*logL + (2*k*n)/(n-k-1);
analysis.BIC  = -2*logL + k*log(n);

params.d1     = x(1);
params.c1     = x(3 : end);
params.s      = x(2);
params.da     = sqrt(2/(1+params.s^2)) * params.s * params.d1;
params.ca     = ( sqrt(2).*params.s ./ sqrt(1+params.s.^2) ) .* params.c1;
params.guess  = guess;

delete current_SDTMLE_data.mat

end


%% log likelihood function
function logL = SDT_logL(parameters)

% loads:
% nR_S1 nR_S2 nCriteria nRatings
load current_SDTMLE_data.mat

% initialize parameters
S1mu = -parameters(1)/2;
S2mu = parameters(1)/2;

S1sd = 1;
S2sd = 1/parameters(2);

c1   = parameters(3 : end);


% calculate response probabilities
ci = [-Inf c1];
cj = [c1 Inf];

pC_S1 = normcdf(cj,S1mu,S1sd) - normcdf(ci,S1mu,S1sd);
pC_S2 = normcdf(cj,S2mu,S2sd) - normcdf(ci,S2mu,S2sd);


pC_S1(pC_S1==0) = 1e-10;
pC_S2(pC_S2==0) = 1e-10;


% calculate likelihood
logL = sum( nR_S1 .* log( pC_S1 ) ) + sum( nR_S2 .* log( pC_S2 ) );

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

logL = -1 * logL;

end