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