function [fits_,sems_,stats_,preds_,resids_] = quickT4_fit(data,guess,lapse)
% function [fits_,sems_,stats_,preds_,resids_] = quickT4_fit(data,guess,lapse)
%
% QUICKT4_FIT fits a weibull function to psychometric data using
%   maximum likelihood maximization. It uses quickT4_err for 
%   error calculation.
%
% 	Input values are:
%     data, in 4 columns...
%		         data(1) = coh  (0 .. 99.9%)
%                data(2) = time (seconds)
%                data(3) = dot dir: left (-1) / right (1)
%                data(4) = correct (1) / error (0)
%
% 	Return values are:
%		fits_   ... see quickT3_val for details
%       	fits_(1) ... A       (coh scale)
%           fits_(2) ... n       (coh exponent)
%           fits_(3) ... B       (time scale)
%       	fits_(4) ... m       (time exponent)
%           fits_(5) ... guess   (guess bias)
%       	fits_(6) ... lapse   ("lapse")
%       sems_   ... Standard errors of the fits. Approximated using the
%                   numerical HESSIAN returned by fmincon (default)
%       stats_  ... [fitLLR Deviance p]
%                   fitLLR is the log likelihood of obtaining the 
%                       data given the fit (returned by quick_err)
%                   Deviance is 2(dataLLR - fitLLR)
%                   p is probability from chi^2 pdf with df = #blocks-3
%       preds_  ... A vector of the probability of making a correct choice given
%                   the fit
%       resids_ ... col 1: log diff between model and actual outcomes
%                   col 2: log diff between model and actual choices (see below)

% Created by jig 9/14/05

if nargin < 1 || isempty(data)
    return
end

global Data

Data = data;

if nargin < 2
    guess = [];
end

if nargin < 3
    lapse = [];
end

%% do the fit
if isempty(guess) && isempty(lapse)
    
    % fit both guess and lapse
    [fits_,f,e,o,l,g,H] = fmincon('quickT4_err', ...
        [ 0.05  0.5 2    0.1   0.0 0.0 ]', [], [], [], [], ...
        [ 0.001 0   0.01 0    -0.4 0.0 ]', ...
        [ 5     3   5    3     0.4 0.5 ]', [], ...
        optimset('LargeScale', 'off', 'Display', 'off', 'Diagnostics', 'off'));
    
elseif isempty(guess)
    
    % fit both guess and lapse
    [fits,f,e,o,l,g,H] = fmincon('quickT4_err', ...
        [ 0.05  0.5 2    0.1  0.0  ]', [], [], [], [], ...
        [ 0.001 0   0.01 0   -0.45 ]', ...
        [ 5     3   5    3    0.45 ]', [], ...
        optimset('LargeScale', 'off', 'Display', 'off', 'Diagnostics', 'off'), ...
        [], lapse);
    fits_  = [fits; lapse];
    
elseif isempty(lapse)
    
    % fit both guess and lapse
    [fits,f,e,o,l,g,H] = fmincon('quickT4_err', ...
        [ 0.05  0.5 2    0.1 0.0 ]', [], [], [], [], ...
        [ 0.001 0   0.01 0   0.0 ]', ...
        [ 5     3   5    3   0.5 ]', [], ...
        optimset('LargeScale', 'off', 'Display', 'off', 'Diagnostics', 'off'), ...
        guess, []);
    fits_  = [fits(1:4); guess; fits(5)];
    
else
    
    % fit both guess and lapse
    [fits,f,e,o,l,g,H] = fmincon('quickT4_err', ...
        [ 0.05  0.5 2    0.1 ]', [], [], [], [], ...
        [ 0.001 0   0.01 0    ]', ...
        [ 5     3   5    3    ]', [], ...
        optimset('LargeScale', 'off', 'Display', 'off', 'Diagnostics', 'off'), ...
        guess, lapse);
    fits_  = [fits; guess; lapse];
end

% Standard errors
%   The covariance matrix is the negative of the inverse of the 
%   hessian of the natural logarithm of the probability of observing 
%   the data set given the optimal parameters.
%   For now, use the numerically-estimated HESSIAN returned by fmincon
%   (which remember is computed from -log likelihood)
if nargout > 1

    % -H because we used -logL in quick_err
    sems_ = sqrt(diag(-((-H)^(-1))));
end

% return stats
if nargout > 2
    
    % log likelihood of the fits ("M1" in Watson)
    % is just the negative of the error function
    M1 = -quickT4_err(fits_);
    
    % deviance is 2(M0 - M1), where
    % M0 is the log likelihood of the data ("saturated model") --
    %   which in this case is, of course, ZERO
    dev = -2*M1;
    
    % probability is from cdf
    if isempty(guess) && isempty(lapse)
        p = 1 - chi2cdf(dev, size(data, 1) - 6);
    elseif isempty(guess) || isempty(lapse)
        p = 1 - chi2cdf(dev, size(data, 1) - 5);
    else
        p = 1 - chi2cdf(dev, size(data, 1) - 4);
    end
    
    stats_ = [M1 dev p];
end

% return the predicted probability (of a correct response) for each observation
if nargout > 3
    preds_ = quickT4_val(Data(:,1:end-1), fits_);
end

% return the deviance residuals ... these are the square roots
%   of each deviance computed individually, signed according to
%   the direction of the arithmatic residual y_i - p_i.
% See Wichmann & Hill, 2001, "The psychometric function: I. Fitting,
%   sampling, and goodness of fit
if nargout > 4

    % resids on correct/error
    cepreds = preds_;
    cepreds(data(:, 4) == 0) = 1 - cepreds(data(:, 4) == 0);
    
    % choice resids are a little trickier - we have to 
    %   make predictions based on choice
    lrpreds = preds_;
    lrpreds(data(:, 3) == -1) = 1 - lrpreds(data(:, 3) == -1);
    choices = data(:, 3);
    choices(data(:, 4) == 0) = -choices(data(:, 4) == 0);
    
    % first column is c/e resids
    % second column is l/r resids
    resids_ = [ ...
        (data(:, 4).*2-1) .* sqrt(-2.*log(cepreds)), ...
        choices .* sqrt(-2.*log(lrpreds))];    
end