function [loglik,m,d] = FitCumGauss_MES(p,r,mset,dset,lapserate,guessrate)
% Fits a cumulative Gaussian to a set of probe values and subject responses
% 
% [loglik,m,d] = FitCumGauss_MES(p,r,mset,dset,lapserate,guessrate)
% Computes likelihood of each combination of PSEs mset and slopes dset
% given a set of probe values and responses using a cumulative Gaussian.
%
% For ease of thinking about it you can view responses as either 0 or 1,
% though in practice anything larger than 0 is treated as 1 and anything
% lower than 0, including 0, is treated as 0. A 1, -1 response scheme as
% input is thus no problem.
%
% By default, a psychometric function ranging from 0% to 100% is used, as
% is suitable for discrimination experiments with a standard in the middle
% of the possible stimulus parameter range. For other paradigms, such as
% detection tasks, one can set the guessrate input to 1/num_alternatives,
% e.g. .5 when doing a 2IFC detection task.
% 
% Optinally fits with a lapse rate, which defaults to 0. If a lapse rate is
% set, the cumulative Gaussian levels off at lapserate/2 and 1-lapserate/2
% when the guessrate is 0. If the guessrate is non zero, the cumulative
% Gaussian that is fit ranges from guessrate to 1-lapserate

% Copyright (c) 2011 by DC Niehorster and JA Saunders

if nargin < 5
    lapserate = 0;
end
if nargin < 6
    guessrate = 0;
end

if guessrate==0
    g0 = lapserate/2;
    g1 = 1 - lapserate;
else
    g0 = guessrate;
    g1 = 1 - lapserate - guessrate;
end
g2 = 1 - g0;

[msamp,dsamp]   = meshgrid(mset,dset);
sz              = size(msamp);
msamp           = msamp(:).';
dsamp           = dsamp(:).';

% the below are equivalent ways of computing the fit, lowest is fastest but
% most obtuse, the others are provided for documentation purposes
if 0
    
    % slow way with a loop,
    % but less memory intensive
    
    loglik = zeros(size(msamp));
    for ksamp = 1:length(msamp)
        for ktrial = 1:length(r)
            if(r(ktrial) > 0)
                currlik = log(g0 + g1*normcdf( (p(ktrial)-msamp(ksamp))/dsamp(ksamp)));
                % reduces to:
                % currlik = log(      normcdf( (p(ktrial)-msamp(ksamp))/dsamp(ksamp)));
                % when lapserate is 0
            else
                currlik = log(g2 - g1*normcdf( (p(ktrial)-msamp(ksamp))/dsamp(ksamp)));
                % currlik = log(1.0 - normcdf( (p(ktrial)-msamp(ksamp))/dsamp(ksamp)));
            end
            loglik(ksamp) = loglik(ksamp) + currlik;
        end
    end
    
elseif 0
    % faster way by blocking, but still looping over trials
    % choose this one if the below get you into lack of memory trouble
    loglik = zeros(size(msamp));
    for ktrial = 1:length(r)
        if(r(ktrial) > 0)
            currlik = log(g0 + g1*normcdf( (p(ktrial)-msamp)./dsamp) );
            % reduces to:
            % currlik = log(      normcdf( (p(ktrial)-msamp)./dsamp));
            % when lapserate is 0
        else
            currlik = log(g2 - g1*normcdf( (p(ktrial)-msamp)./dsamp) );
            % currlik = log(1.0 - normcdf( (p(ktrial)-msamp)./dsamp));
        end
        loglik     = loglik + currlik;
    end
    
elseif 0
    % faster way by blocking,
    % but more memory intensive
    rmat        = repmat(r(:)>0,[1,length(msamp)]);
    mmat        = repmat(msamp,[length(r),1]);
    dmat        = repmat(dsamp,[length(r),1]);
    pmat        = repmat(p(:),[1,length(msamp)]);
    
    temp        = g2 - g1*normcdf( (pmat-mmat)./dmat);
    ind         = find(rmat(:));
    temp(ind)   = 1 - temp(ind);
    loglik      = sum(log(temp),1);
    
else
    % slightly faster and less memory intensive way than previous, using
    % bsxfun (very slightly actually, but hey)
    temp        = g2 - g1*normcdf( bsxfun(@rdivide,bsxfun(@minus,p(:),msamp),dsamp));
    rmat        = repmat(r(:)>0,[1,length(msamp)]);
    temp(rmat)  = 1 - temp(rmat);
    loglik      = sum(log(temp),1);
end

[~,kmax]= max(loglik);
m       = mean(msamp(kmax));
d       = mean(dsamp(kmax));

loglik  = reshape(loglik,sz);

%imagesc(exp(0.5*loglik))