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

% fit = fit_rs_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 RESPONSE-SPECIFIC type 2 SDT analysis of the data. i.e.
% meta-d' is computed separately for each response type. 
%
% N.B. it is still important to provide input including data from both 
% response types, even if you are only interested in the analysis for one 
% of the response types.
% 
% 
%
% 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
%
% N.B. if nR_S1 or nR_S2 contain zeros, this may interfere with estimation of
% meta-d'.
%
% Some options for dealing with response cell counts containing zeros are:
% 
% (1) Add a small adjustment factor, e.g. adj_f = 1/(length(nR_S1), to each 
% input vector:
% 
% adj_f = 1/length(nR_S1);
% nR_S1_adj = nR_S1 + adj_f;
% nR_S2_adj = nR_S2 + adj_f;
% 
% This is a generalization of the correction for similar estimation issues of
% type 1 d' as recommended in
% 
% Hautus, M. J. (1995). Corrections for extreme proportions and their biasing 
%     effects on estimated values of d'. Behavior Research Methods, Instruments, 
%     & Computers, 27, 46-51.
%     
% When using this correction method, it is recommended to add the adjustment 
% factor to ALL data for all subjects, even for those subjects whose data is 
% not in need of such correction, in order to avoid biases in the analysis 
% (cf Snodgrass & Corwin, 1988).
% 
% (2) Collapse across rating categories.
% 
% e.g. if your data set has 4 possible confidence ratings such that length(nR_S1)==8,
% defining new input vectors
% 
% nR_S1_new = [sum(nR_S1(1:2)), sum(nR_S1(3:4)), sum(nR_S1(5:6)), sum(nR_S1(7:8))];
% nR_S2_new = [sum(nR_S2(1:2)), sum(nR_S2(3:4)), sum(nR_S2(5:6)), sum(nR_S2(7:8))];
% 
% might be sufficient to eliminate zeros from the input without using an adjustment.
%
% * 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.
% For most purposes, we recommend setting s = 1. 
% See http://www.columbia.edu/~bsm2105/type2sdt for further discussion.
%
% * 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.t1ca      = type 1 criterion for overall data, RMS units
% fit.s         = sd(S1) / sd(S2)
%
% fit.meta_da_rS1   = meta-d' for S1 responses, RMS units
% fit.t1ca_rS1      = type 1 criteron for meta-d' fit for S1 responses, RMS units
% fit.t2ca_rS1      = type 2 criteria for meta-d' fit for S1 responses, RMS units
% fit.M_ratio_rS1   = meta_da_rS1 / da
% fit.M_diff_rS1    = meta_da_rS1 - da
%
% fit.logL_rS1      = log likelihood of meta-d' fit for S1 responses
%
% fit.obs_HR2_rS1   = actual type 2 hit rates for S1 responses
% fit.est_HR2_rS1   = estimated type 2 hit rates for S1 responses
% fit.obs_FAR2_rS1  = actual type 2 false alarm rates for S1 responses
% fit.est_FAR2_rS1  = estimated type 2 false alarm rates for S1 responses
%
%
% fit.meta_da_rS2   = meta-d' for S2 responses, RMS units
% fit.t1ca_rS2      = type 1 criteron for meta-d' fit for S2 responses, RMS units
% fit.t2ca_rS2      = type 2 criteria for meta-d' fit for S2 responses, RMS units
% fit.M_ratio_rS2   = meta_da_rS2 / da
% fit.M_diff_rS2    = meta_da_rS2 - da
%
% fit.logL_rS2      = log likelihood of meta-d' fit for S2 responses
%
% fit.obs_HR2_rS2   = actual type 2 hit rates for S2 responses
% fit.est_HR2_rS2   = estimated type 2 hit rates for S2 responses
% fit.obs_FAR2_rS2  = actual type 2 false alarm rates for S2 responses
% fit.est_FAR2_rS2  = estimated type 2 false alarm rates for S2 responses
%
%
% 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.

% 2015/07/23 - added comments to help section as well as a warning output 
%              for nR_S1 or nR_S2 inputs containing zeros

%% parse inputs

% check 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 any(nR_S1 == 0) || any(nR_S2 == 0)
    disp(' ')
    disp('WARNING!!')
    disp('---------')
    disp('Your inputs')
    disp(' ')
    disp(['nR_S1 = [' num2str(nR_S1) ']'])
    disp(['nR_S2 = [' num2str(nR_S2) ']'])
    disp(' ')
    disp('contain zeros! This may interfere with proper estimation of meta-d''.')
    disp('See ''help fit_rs_meta_d_MLE'' for more information.')
    disp(' ')
    disp(' ')
end

% assign input default values
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;


%% find actual type 2 FAR and HR (data to be fit)

% 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
    t2FAR_rS2(i-1) = sum( I_nR_rS2(i:end) ) / sum(I_nR_rS2);
    t2HR_rS2(i-1)  = sum( C_nR_rS2(i:end) ) / sum(C_nR_rS2);
    
    t2FAR_rS1(i-1) = sum( I_nR_rS1(i:end) ) / sum(I_nR_rS1);
    t2HR_rS1(i-1)  = sum( C_nR_rS1(i:end) ) / sum(C_nR_rS1);    
end

t2FAR_rS1 = t2FAR_rS1(end:-1:1);
t2HR_rS1  = t2HR_rS1(end:-1:1);



%% set up constraints

% parameters
% meta-d' - 1
% t2c     - (nCriteria-1)/2

A = [];
b = [];

nCritPerResp = (nCriteria-1)/2;

offset = 1;
for crit = 1 : nCritPerResp - 1
    % c(crit) <= c(crit+1) --> c(crit) - c(crit+1) <= .001
    A(end+1,[offset+crit offset+crit+1]) = [1 -1];
    b(end+1) = -.001;
end

LB_rS1 = [];
LB_rS1 = [LB_rS1 -10];                      % meta-d' rS1
LB_rS1 = [LB_rS1 -20*ones(1,nCritPerResp)]; % criteria lower than t1c

UB_rS1 = [];
UB_rS1 = [UB_rS1 10];                       % meta-d' rS1
UB_rS1 = [UB_rS1 zeros(1,nCritPerResp)];    % criteria lower than t1c



LB_rS2 = [];
LB_rS2 = [LB_rS2 -10];                      % meta-d' rS2
LB_rS2 = [LB_rS2 zeros(1,nCritPerResp)];    % criteria higher than t1c

UB_rS2 = [];
UB_rS2 = [UB_rS2 10];                       % meta-d' rS2
UB_rS2 = [UB_rS2 20*ones(1,nCritPerResp)];  % criteria higher than t1c



%% select constant criterion type

constant_criterion_rS1 = 'meta_d1_rS1 * (t1c1 / d1)'; % relative criterion
constant_criterion_rS2 = 'meta_d1_rS2 * (t1c1 / d1)'; % relative criterion


%% set up guess

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_rS1 = d1;
meta_d1_rS2 = d1;

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


guess_rS1 = [meta_d1_rS1 t2c1(1:nCritPerResp) - eval(constant_criterion_rS1)];
guess_rS2 = [meta_d1_rS2 t2c1(nCritPerResp+1 : end) - eval(constant_criterion_rS1)];

save fitM_data.mat nR_S1 nR_S2 t2FAR_rS2 t2HR_rS2 t2FAR_rS1 t2HR_rS1 nRatings t1c1 s d1 fncdf constant_criterion_rS1 constant_criterion_rS2


%% fit for S1 responses

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

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

[x f]=fmincon(@fitM_rS1_logL,guess_rS1,A,b,[],[],LB_rS1,UB_rS1,[],op);

meta_d1_rS1 = x(1);
meta_c1_rS1 = eval(constant_criterion_rS1);

t2c1_rS1    = x(2:end) + eval(constant_criterion_rS1);

logL_rS1    = -f;


%% find model-estimated t2FAR and t2HR for S1 responses

% find the estimated type 2 FAR and HR
S1mu = -meta_d1_rS1/2; S1sd = 1;
S2mu =  meta_d1_rS1/2; S2sd = S1sd/s;

% adjust so that everything is centered on t1c1=0
h = 1-normcdf(0,S2mu,S2sd);
f = 1-normcdf(0,S1mu,S1sd);

shift_c1 = (-1/(1+s)) * (norminv(h)+norminv(f)); % this is the value of c1 midway b/t S1 and S2

S1mu = S1mu + shift_c1; % shift S1 and S2mu so that they lie on an axis for 0 --> c1=0
S2mu = S2mu + shift_c1;


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


for i=1:length(t2c1_rS1)
    
    t2c1_lower = t2c1_rS1(i);
        
    I_FAR_area_rS1 = fncdf(t2c1_lower,S2mu,S2sd);
    C_HR_area_rS1  = fncdf(t2c1_lower,S1mu,S1sd);
    
    est_t2FAR_rS1(i) = I_FAR_area_rS1 / I_area_rS1;
    est_t2HR_rS1(i)  = C_HR_area_rS1 / C_area_rS1;
    
end



%% fit for S2 responses

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

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

[x f]=fmincon(@fitM_rS2_logL,guess_rS2,A,b,[],[],LB_rS2,UB_rS2,[],op);

meta_d1_rS2 = x(1);
meta_c1_rS2 = eval(constant_criterion_rS2);

t2c1_rS2    = x(2:end) + eval(constant_criterion_rS2);

logL_rS2    = -f;


%% find estimated t2FAR and t2HR

% find the estimated type 2 FAR and HR
S1mu = -meta_d1_rS2/2; S1sd = 1;
S2mu =  meta_d1_rS2/2; S2sd = S1sd/s;

% adjust so that everything is centered on t1c1=0
h = 1-normcdf(0,S2mu,S2sd);
f = 1-normcdf(0,S1mu,S1sd);

shift_c1 = (-1/(1+s)) * (norminv(h)+norminv(f)); % this is the value of c1 midway b/t S1 and S2

S1mu = S1mu + shift_c1; % shift S1 and S2mu so that they lie on an axis for 0 --> c1=0
S2mu = S2mu + shift_c1;


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


for i=1:length(t2c1_rS2)
    
    t2c1_upper = t2c1_rS2(i);
        
    I_FAR_area_rS2 = 1-fncdf(t2c1_upper,S1mu,S1sd);
    C_HR_area_rS2  = 1-fncdf(t2c1_upper,S2mu,S2sd);
       
    est_t2FAR_rS2(i) = I_FAR_area_rS2 / I_area_rS2;
    est_t2HR_rS2(i)  = C_HR_area_rS2 / C_area_rS2;
    
end


%% package output

% type 1 params
fit.da         = SDT_s_convert(d1,   s,'d1','da');
fit.t1ca       = SDT_s_convert(t1c1, s,'c1','ca');
fit.s          = s;


% type 2 fits for rS1
fit.meta_da_rS1   = SDT_s_convert(meta_d1_rS1, s,'d1','da');
fit.t1ca_rS1      = SDT_s_convert(meta_c1_rS1, s,'c1','ca');
fit.t2ca_rS1      = SDT_s_convert(t2c1_rS1,    s,'c1','ca');

fit.M_ratio_rS1   = fit.meta_da_rS1 / fit.da;
fit.M_diff_rS1    = fit.meta_da_rS1 - fit.da;

fit.logL_rS1      = logL_rS1;

fit.obs_HR2_rS1   = t2HR_rS1;
fit.est_HR2_rS1   = est_t2HR_rS1;
fit.obs_FAR2_rS1  = t2FAR_rS1;
fit.est_FAR2_rS1  = est_t2FAR_rS1;


% type 2 fits for rS2
fit.meta_da_rS2   = SDT_s_convert(meta_d1_rS2, s,'d1','da');
fit.t1ca_rS2      = SDT_s_convert(meta_c1_rS2, s,'c1','ca');
fit.t2ca_rS2      = SDT_s_convert(t2c1_rS2,    s,'c1','ca');

fit.M_ratio_rS2   = fit.meta_da_rS2 / fit.da;
fit.M_diff_rS2    = fit.meta_da_rS2 - fit.da;

fit.logL_rS2      = logL_rS2;

fit.obs_HR2_rS2   = t2HR_rS2;
fit.est_HR2_rS2   = est_t2HR_rS2;
fit.obs_FAR2_rS2  = t2FAR_rS2;
fit.est_FAR2_rS2  = est_t2FAR_rS2;


% S1 units
fit.S1units.d1   = d1;
fit.S1units.t1c1 = t1c1;

fit.S1units.meta_d1_rS1 = meta_d1_rS1;
fit.S1units.t1c1_rS1    = meta_c1_rS1;
fit.S1units.t2c1_rS1    = t2c1_rS1;

fit.S1units.meta_d1_rS2 = meta_d1_rS2;
fit.S1units.t1c1_rS2    = meta_c1_rS2;
fit.S1units.t2c1_rS2    = t2c1_rS2;

delete fitM_data.mat

end


%% function to find the SSE for type 2 HR and FAR for S1 responses
function logL = fitM_rS1_logL(input)

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

% load the behavioral data
% t2FAR_rS2 t2HR_rS2 t2FAR_rS1 t2HR_rS1 nRatings t1c1 s d1 fncdf constant_criterion_rS1 constant_criterion_rS2
load fitM_data.mat


% find the estimated type 2 FAR and HR
S1mu = -meta_d1_rS1/2; S1sd = 1;
S2mu =  meta_d1_rS1/2; S2sd = S1sd/s;


% adjust so that everything is centered on t1c1=0
h = 1-normcdf(0,S2mu,S2sd);
f = 1-normcdf(0,S1mu,S1sd);

shift_c1 = (-1/(1+s)) * (norminv(h)+norminv(f)); % this is the value of c1 midway b/t S1 and S2

S1mu = S1mu + shift_c1; % shift S1 and S2mu so that they lie on an axis for 0 --> c1=0
S2mu = S2mu + shift_c1;


% adjust so that t1c1 = 0
S1mu = S1mu - eval(constant_criterion_rS1);
S2mu = S2mu - eval(constant_criterion_rS1);
t1c1 = 0;



%%% set up MLE analysis

% get type 2 response counts

for ic = 1:nRatings
    
    % S1 responses
    nC_rS1(ic) = nR_S1(ic);
    nI_rS1(ic) = nR_S2(ic);
        
end

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


t2c1x = [-Inf t2c1 t1c1];

% for ic = 1:4
for ic = 1:nRatings
    prC_rS1(ic) = ( fncdf(t2c1x(ic+1),S1mu,S1sd) - fncdf(t2c1x(ic),S1mu,S1sd) ) / C_area_rS1;
    prI_rS1(ic) = ( fncdf(t2c1x(ic+1),S2mu,S2sd) - fncdf(t2c1x(ic),S2mu,S2sd) ) / I_area_rS1;   
end
    

logL = 0;
for ic = 1:nRatings
    logL = logL + nC_rS1(ic)*log(prC_rS1(ic)) + nI_rS1(ic)*log(prI_rS1(ic));
end

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

logL = -logL;

end


%% function to find the SSE for type 2 HR and FAR for S2 responses
function logL = fitM_rS2_logL(input)

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

% load the behavioral data
load fitM_data.mat


% find the estimated type 2 FAR and HR
S1mu = -meta_d1_rS2/2; S1sd = 1;
S2mu =  meta_d1_rS2/2; S2sd = S1sd/s;


% adjust so that everything is centered on t1c1=0
h = 1-normcdf(0,S2mu,S2sd);
f = 1-normcdf(0,S1mu,S1sd);

shift_c1 = (-1/(1+s)) * (norminv(h)+norminv(f)); % this is the value of c1 midway b/t S1 and S2

S1mu = S1mu + shift_c1; % shift S1 and S2mu so that they lie on an axis for 0 --> c1=0
S2mu = S2mu + shift_c1;


% adjust so that t1c1 = 0
S1mu = S1mu - eval(constant_criterion_rS2);
S2mu = S2mu - eval(constant_criterion_rS2);
t1c1 = 0;



%%% set up MLE analysis

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


% get type 2 probabilities
C_area_rS2 = 1-fncdf(t1c1,S2mu,S2sd);
I_area_rS2 = 1-fncdf(t1c1,S1mu,S1sd);

t2c1x = [t1c1 t2c1 Inf];

for ic = 1:nRatings
    prC_rS2(ic) = ( (1-fncdf(t2c1x(ic),S2mu,S2sd)) - (1-fncdf(t2c1x(ic+1),S2mu,S2sd)) ) / C_area_rS2;
    prI_rS2(ic) = ( (1-fncdf(t2c1x(ic),S1mu,S1sd)) - (1-fncdf(t2c1x(ic+1),S1mu,S1sd)) ) / I_area_rS2;
end
    

logL = 0;
for ic = 1:nRatings
    logL = logL + nC_rS2(ic)*log(prC_rS2(ic)) + nI_rS2(ic)*log(prI_rS2(ic));
end

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

logL = -logL;


end



function out = SDT_s_convert(in1,s,cont1,cont2)

% out = SDT_s_convert(in1,s,cont1,cont2)
%
% in1 - input var to be converted
% s - sd(S1)/sd(S2)
% cont1 - identity of in1
% cont2 - identity of out
%
% cont1 and cont2 can be these tokens
% 'da','d1','d2','ca','c1',c2'

switch cont1
    
    % convert d'
    % s = d2 / d1
    % da = sqrt(2./(1+s.^2)) .* d2
    case 'da'
        da = in1;
        d2 = da ./ sqrt(2./(1+s.^2));
        d1 = d2 ./ s;
        
    case 'd1'
        d1 = in1;
        d2 = d1 .* s;
        da = sqrt(2./(1+s.^2)) .* d2;
        
    case 'd2'
        d2 = in1;
        d1 = d2 ./ s;
        da = sqrt(2./(1+s.^2)) .* d2;
        
    % convert c
    % s = c2 / c1
    % ca = ( sqrt(2).*s ./ sqrt(1+s.^2) ) .* c1;
    case 'ca'
        ca = in1;
        c1 = ( sqrt(1+s.^2) ./ sqrt(2).*s ) .* ca;
        c2 = c1 .* s;
        
    case 'c1'
        c1 = in1;
        c2 = c1 .* s;
        ca = ( sqrt(2).*s ./ sqrt(1+s.^2) ) .* c1;
        
    case 'c2'
        c2 = in1;
        c1 = c2 ./ s;
        ca = ( sqrt(2).*s ./ sqrt(1+s.^2) ) .* c1;
end

eval(['out = ' cont2 ';']);

end