% function [bestparams, statstruct ] =
% psyfitnew(datasource, nchoice, <exponent/mode, default 0> , <interactive, default 1> )
%
% ***********ARGUMENTS: one argument, datasource, is required.
%
% DATASOURCE can be a file name string, or a matrix in the workspace. If
% a file, it should be an ascii file with just the data matrix. Data matrix
% must have 3 or 4 columns:
% Col 1: a list of tested stimulus values.
% Col 2: number of trials presented (at that value for this row of data).
% Col 3: number of positive responses "yes/correct" in those trials.
% Col 4 (optional) a condition index: must be 1 thru N for N conditions. 
%
% Col 1 tested values can be different for each condition, and 
% Col 1 values can be repeated indefinitely in any order, e.g. if
% each data row represents a single trial (so Col 2 entry = 1)
% or if multiple data files are merged by concatenating their data. 
%
%
% OPTIONAL ARGUMENTS:
%
% NCHOICE 
% nchoice = 1: for YES/NO (default); 
%   fits full normal ogive, 0 < p < 1 for stimval = -Inf:Inf;
%   Ceiling on p is 1-LAPSEPROB, where
%   LAPSEPROB is set in the script, currently 0.02.
% nchoice >= 2: for nAFC FORCED-CHOICE; floor for p is 1/n.
% nchoice = 0: for COMPARATIVE JUDGMENT, fits full normal ogive.
%   To accommodate finger errors, fitted p ranges from  
%    fingerrprob/2 to (1 - fingerrprob/2) if nchoice =0. 
%    fingerrprob is set in the script, currently 0.02.
%   Otherwise equivalent to nchoice = 1.
%
% EXPONENT/MODE, if supplied and non-zero, requests POWER_OF_X mode.
% Normally, a full normal ogive is fit to the proportion correct as f(X),
% with chance performance at X = -Inf. In POWER_OF_X mode on the other
% hand, only the top half of the normal cdf is used for the fit, and the
% inflection usually required is introduced by transforming to X^exponent.
%
% In default mode (exponent/mode not supplied or = 0), the stimulus values
% X could be negative (e.g. logs of some magnitude). But in POWER_OF_X
% mode, the exponent is the power n of the stimulus value X, where pcorrect
% is the top half of a normal cdf as a function of X^n, and X must be all
% positive. POWER_OF_X mode is used if and only if exponent is nonzero.
% This has a theoretical rationale for 2AFC where it corresponds to
% predictions of a model where Gaussian noise is added to a signal that
% increases as X^exponent. For other cases, the fitted curve is pinned at X
% = 0 to chance = n/nchoice, or to 0 +ve responses if nchoice = 1 (for
% yes/no or free response experiments). The second parameter is always the
% sigma of the associated normal cdf, that is the step in stimulus value
% from 50% to 84% correct as nchoice _.Inf, from (.5 + .5/nchoice) to (.84
% +.84/nchoice) in the general case.
% If the exponent argument is negative, it makes no sense to adopt it, 
% and the exponent is iteratively optimized instead.
%
% The exponent choice determines the steepness of the fitted curve for a
% given threshold value, and threshold in this mode is the only free
% parameter if the exponent supplied as an argument is adopted. But the
% code does allow for the exponent to be varied as a second parameter.
% Presently this option is only available by user choice if the final
% argument, interactive, is non-zero.
%
% MAKEPLOTS = 1 (default)for plotting final results, showing control
% options; 0 for no plots; > 1 to plot intermediate results as well.
% Psyfitnew drops two of the esoteric plots available in psyfit.
%
% ***********USAGE:
% [bestparams, statstructarray] = psyfit(datasource)
% fits a full normal ogive, 0 < p < 1.
% [bestparams, statstructarray] = psyfit(datasource,4)
% fits a full normal ogive with a chance (25% correct) asymptote at
% –Infinity stimulus value. This is a natural choice when (and only when)
% the stimulus is the log of a magnitude, since above-chance performance is
% predicted at zero stimulus value.
% 
% POWER_OF_X mode requires a non-zero exponent/mode argument:
% [bestparams, statstructarray] = psyfit(datasource,4,3,0)
% or
% [bestparams, statstructarray] = psyfit(datasource,4,3)
% will use POWER_OF_X mode, fixing the exponent at 3 by default with no 
% opportunity for the user to free the exponent parameter.
% But
% [bestparams, statstructarray] = psyfit(datasource,4,3,1)
% allows the user to abandon the default exponent value so it becomes a
% free second parameter for best fit. This is helpful since the syntax 
% requires a nonzero value to trigger this mode, but the exponent
% value may not be one that you have confidence in.
%
%
% *********************OUTPUTS:
%
% BESTPARAMS reports the fit results as (generally) 2 parameter values.
% 
% For yes/no or comparative judgments (nchoice=1 or 0), the first parameter
% is the 50% point or inflection point of the fitted ogive. This is the
% standard definition of threshold in yes/no experiments and serves equally
% well as an estimate of the point of subjective equality in a comparative
% judgment (or matching) experiment. With nAFC, param is where p is halfway
% between chance and max. The second parameter is the standard
% deviation of the cumulative normal on which the fit is based.
% 
% For forced-choice data analysed in POWER_OF_X mode, the first parameter
% is the standard deviation of the cumulative normal on which the fit is
% based, and the second parameter is the exponent on stimulus value (if
% freely varied rather than fixed as argument 3) as described above.
% 
% The output 'bestparams' is a row vector if there is only one condition.
% Generally it is a matrix with one row per condition. For a two parameter
% fit it has four columns: columns 1 and 3 are max likelihood threshold and
% its standard error. Columns 2 and 4 are the sigma or exponent
% (in POWER_OF_X mode) parameter and its standard error. For a one parameter fit
% (i.e. the exponent argument in POWEROF_X mode is supplied) there are
% only the first 2 columns: mean threshold and its standard error.
% 
% STATSTRUCT 
% If there are multiple conditions (datamatrix has 4 cols), the output
% 'statstruct {ic}' is an structure array, one structure per condition,
% with fields for results, including bestparams; fields of 'statstruct {ic}'
% are generally scalar or 2 long rows with one entry per free parameter.
% statstruct{ic}.marglh and statstruct {ic}.pvals define the 2  marginal
% pdfs, 2 by N. The 'fits' field of statstruct{ic} is a matrix of size
% nxvals(ic) by 4, where nxvals is the number of xvalues actually presented
% for that condition. The 4 columns show the xvalue, number of
% presentations, observed proportion correct and predicted proportion
% correct. Thus statstruct {ic}.fits =[xvalues ntotal observedpcorrect
% predictedpcorrect].
%
% Various statistics are given...see 'ESTIMATES' below. Notably, meanlatticemarg
% is the mean of the marginal pdf for the value of parameter 1 (column 1)
% and parameter 2 (column 2). sigmalatticemarg is the s.d of each pdf.
% Quantiles are given at quantileLevels = [.005 .025 .05 .5 .95 .975 .995],
% for threshold (row 1) and for slope parameter (row 2). These define
% 1%, 5% and 10% confidence intervals plus the median of the pdfs.
% The marglh element of statstruct is the probability of the data for the
% sampled values of each  model parameter, using empirical Bayes weighting
% over the free parameter. The sampled values are platticevals in the same
% 2 by nvals format as marglh, and quantiles are derived by interpolation
% in marglh. In statstruct, the bestparams element does not contain the 
% standard errors that are included in the separately returned bestparams
% output argument. Instead it contains only two scalars: the mean and sigma
% of the MAP fitted normal cdf.
%
% If there is only one condition (datamatrix has only 3 cols)
% statstruct is a single structure rather than the array above.
%
%
% ****************ALGORITHM:
% Fits a cumulative Gaussian to data from binary response experiments in
% Nxvals x 3 or 4 data format. Two basic modes, selected by choice of
% nchoice argument:
% 1a) FULL OGIVE, YES/NO OPTION (nchoice = 1) Gaussian mean
% bestparams(1) and std dev (bestparams(2))are varied. Here x is often
% a log value, with response probability zero at -Inf. THe returned
% threshold is the 'half-way up' point on the ogive, i.e. 50% correct
% 1b) FULL OGIVE, nAFC
% option (nchoice = n) fits as in 1 but scales the predicted proportion
% correct to a range from 1/nchoice to 1 (allowing for lapses in the usual
% way). Again x is typically a LOG value, giving chance performance
% (proportion correct = 1/nchoice) at x = -Inf. Again the returned
% threshold is the 'half-way up' point on the ogive, i.e. 50% of the way
% from 1/nchoice correct to 100% correct.
% 2) POWER_OF_X OPTION
% With this option, usually invoked for nAFC, a HALF normal ogive is
% fitted, with a scaling threshold x0 (bestparams(1)) and n
% (bestparams(2))a best fitting exponent on x . Thus predicted proportion
% correct is normcdf(x/x0)^n. In this case x has to be LINEAR POSITIVE
% intensity and chance performance is expected at x = 0. The threshold x0
% (= bestparams(1)) corresponds to a standard normal deviate of +1, or to
% 84% correct. In this mode, the exponent can be fixed at the value of the
% third, 'exponent' argument...useful since the exponent is hard to
% estimate from limited data. The implicit model is a POWER_OF_X threshold
% nonlinearity in the response as a function of x, followed by Gaussian
% noise added to the response (Kontsevich and Tyler).
% For y/n, the half ogive is doubled in range to reach chance=.5 at X = 0.
% For FC with nchoice not 2, the range is set to reach chance (1/nchoice)
% at X=0.
% The use of POWER_OF_X mode with ragged data
% and with variable exponent may be troublesome. In that mode, it is safest
% to fix the exponent by supplying it as an argument. As a compromise,
% a prior for slope can be applied in POWER_OF_X mode within the halfCDFllh
% function. This prior is centered at a default exponent of 2.5, with sigma
% around 40%, see halfCDFllh to change.
%
% Sorting by xval is done first, then maximum likelihood fits are done to
% the frequency-of-seeing data, by calling fitrt (fullCDFllh.m or
% halfCDFllh.m) to return minus the log likelihood, and to append log 
% likelihood and associated parameter values to the list LOGLHS for each 
% iteration. 
%
% **************ESTIMATES:
% 1. The maximum likelihood point (2D mode of the likelihood distribution,
% by quadratic interpolation) is internally returned by fminsearch in
% optima. For a good data set, this is a good estimate, but (1) the
% center of gravity of the likelihood distribution has some advantages,
% especially that it integrates over the slope parameter so as to take
% account of the uncertainty in that in the estimate of threshold. This
% center of gravity is obtained by combining the mean of the marginal
% distributions. The mean minimizes the expected loss if errors in
% parameter estimation have a cost proportional to their squares. This is
% one reason to prefer it to the mode. A direct mean and s.d. of the
% marginal lattice pdf values is not a good idea though, since these are
% heavily influenced by truncation and sampling, and are completely
% inappropriate if sampling is not uniform. So (2) a quadratic fit to the
% near-peak marginals avoids this problem and may be preferable. Such a fit is
% what generates meanquad and sigmaquad and those estimators are selected if
% USEQUAD is nonzero in line 475 or so. But a problem with the quadratic
% fit is that confidence intervals should be constrained by (only) the
% flanks of the pdf, and the peak only by the highest points (of which only
% 3 are needed for interpolation with a Gaussian approximation). Therefore
% (3) Bayes; doesn't indicate a proper weighting, and
% suggests instead getting CIs rather than sigmas by interpolation in the
% marginal pdf. The CIS are crude, from cumtrapz; could use a denser array of
% marginal pdf values, obtained by cubic interpolation in the log of the
% margpdf, much as in FastEstimate, but with prior interpolation so that
% the intervals spanned by the trapezoids are less than for the initial
% lattice.
%
% Finally the prelimsigmas are used only to set up spacing of the test 
% lattice in parameter space. They are based on 1D parameter variation. 
% They assume independence of the 2D pdf in parameter space, which
% can sometimes help in handling difficult cases, but parameter dependence
% calls for some inflation of the prelim estimates, by the empirically
% chosen (and non-critical) constant PRELIMSIGMASCALER. 
% The prelimsigmas are based the 1D cross-sections of likelihood
% distribution through the 2D optimum in parameter space. The peak 
% curvature and SD of the cross-section are each evaluated, and the larger
% value is adopted, for reasons stated in notes below.  
%

% don macleod (mainly: some bits based on ed vul's FAST code).
%****************************************************************************

function [bestparams, statstruct] = psyfit(datasource, nchoice, exponent, interactive )
clear global; 
global ORTHOG;
global TWOPARAM DEFAULTEXP NCHOICE;
global XVALS NTT NC LOGLHS ITERATION;
%% Ask user for data source and relevant column indices as needed:
MAKEPLOTS = 1; % default; set >1 for extra plots, zero for none.
if nargin==0 % have to ask for data file as source
    % File assumed to be in current directory; consider uigetfile()
    datasource = input('ascii file to analyze, including any extension: ','s');
elseif nargin < 2
    NCHOICE = 1;  % default; or ask: input('Enter 1 if experiment was yes/no detection, Nchoices if it was forced choice, or 0 for comparative judgment:  ');
elseif nargin >= 2
    NCHOICE = nchoice;
    if nargin >= 4 % overwrite default for output control
        MAKEPLOTS = interactive;
    end
end


TWOPARAM='Y'; % set to 'Y' to obtain 2 parameter f/c fits in 2AFC POWER_OF_X mode;'N' for one parameter
% (varies threshold + exponent if 'Y', threshold only if 'N'.)
PRELIMSIGMASCALER = 2; % guessed compensation for parameter dependence
fitrt = 'fullCDFllh'; % Usual mode, use fullCDFllh below fit full ogive to x ranging from -Inf to +Inf
if nargin >= 3  % POWER_OF_X mode is wanted if a non-zero exponent has been supplied...
    if exponent~=0 
        if nchoice<2
            warning('error: nonzero exponent (POWER_OF_X mode) with nchoice < 2');  % better: error?
        end
        DEFAULTEXP = exponent;  % use the exponent passed as an argument...
        fitrt='halfCDFllh'; % POWER_OF_X mode by calling halfCDFllh below
        TWOPARAM = 'N';  % use the passed exponent argument unless chosen otherwise interactively below
        if exponent>0 && MAKEPLOTS>=1 % if interactive, check default exponent is OK
            disp('Using 2AFC POWER_OF_X mode for linear positive stimulus values, threshold will be = 1 s.d. of normal cdf, or 84% correct');
            exponent % if interactive, allow user to opt for optimization of the exponent
            tmpstring = input('POWER_OF_X mode: accept above value (y) or vary exponent for best fit (n)?  ', 's');
            if (tmpstring(1)=='n' || tmpstring(1)=='N')
                TWOPARAM = 'Y';  
            end
        end
    end
    if exponent < 0
        TWOPARAM = 'Y'; % iteratively optimize, never use DEFAULTEXP
    end
end


%% ************  begin with N x 3 or 4 matrix, in workspace or from ascii data file:
if ~ischar(datasource) % workspace
    allthedata = datasource;
else
    allthedata = load(datasource); % file
end
[h w]=size(allthedata);
if all(fitrt=='halfCDFllh') && any(datasource(:,1) < 0)
    error('POWER_OF_X mode not available when stimulus values can be negative; set exponent to zero, or omit it, for a regular full ogive fit');
end

if w > 3
    condindices = unique(allthedata(:,4));
    allthedata = sortrows(allthedata, [1 4]);
else
    condindices = 1;
    allthedata = sortrows(allthedata, 1);
    allthedata(:,4) = ones(size(allthedata(:,1)));  % fill 4th column with ones...simpler later
end
nconds=length(condindices);

%%  Consolidate the data for stimuli presented in each condition
for ic = 1:nconds
    % Put NaNs in statstruct{ic} structure as defaults in cases not defined below
    statstruct{ic}.bestparams = NaN.*ones(1,2);
    statstruct{ic}.prelimsigmas = NaN.*ones(1,2);  % preliminary estimates, superseded by sigmaquad etc
    statstruct{ic}.meanlatticemarg = NaN.*ones(1,2); % mean of raw lattice sample points from marginal pdf(s)
    statstruct{ic}.sigmalatticemarg = NaN.*ones(1,2);
    statstruct{ic}.meanquad = NaN.*ones(1,2);  % mean of gaussian fit(s) to 1d marginal pdf(s)
    statstruct{ic}.sigmaquad = NaN.*ones(1,2); % sd of gaussian(s)
    disp(['Condition ' num2str(ic)]);
    ITERATION=1;
    data_thiscond = allthedata(allthedata(:,4) == ic,:); % relevant rows, even if no presentations
    data_thiscond = data_thiscond(data_thiscond(:,2) >= 0,:); % rows for presented stimulus values
    XVALS = sort(unique(data_thiscond(:,1)));
    for xindex = 1:length(XVALS)
        NTT(xindex,1) = sum(data_thiscond(data_thiscond(:,1)==XVALS(xindex,1),2)); %ntotal
        NC(xindex,1) =  sum(data_thiscond(data_thiscond(:,1)==XVALS(xindex,1),3)); % %NCorrect
    end
    k = find(NTT(1:xindex) > 0); % restrict to range of x presented in current condition
    NTT = NTT(k);
    XVALS = XVALS(k);
    NC = NC(k);
    Npresented{ic} = NTT;
    xpresented{ic} = XVALS;
    opc{ic} = NC./NTT; %observed prop correct for condition, indexed by ic; for plot vs xpresented
    xused = xpresented{ic}(isfinite(xpresented{ic})) ; % & opc{ic}>0 & opc{ic} < 1); % truncate to define sigma of stimuli presented with 0 < p < 1?
    sigmastimvals{ic} = sqrt(var( xused, opc{ic}(isfinite(xpresented{ic})) )); % s.d. of x presented; refine? /2? 
%% Initial iterations for max likelihood estimates of the parameters:call fitrt to get those
    % Note:fitrt appends to LOGLHS = list of llhs with associated paramvals
    meanval=sum(XVALS.*NTT)./sum(NTT);
    disp('Iterating...')
    if all(fitrt=='halfCDFllh'),
        initparams=[meanval log(2)]; % guessed scaler and log_e exponent for POWEROFX model; don't choose log(1)...
        if TWOPARAM=='N'
            initparams=meanval;
        end
        optima = fminsearch(@halfCDFllh, initparams); %  warns, keeps initparams as initial guess if maxiter exceeded
        if max(LOGLHS(:,1)-LOGLHS(end,1) > 0.1)
            [tmp, tmpi] = max(LOGLHS(:,1));
            optima = fminsearch(@halfCDFllh, LOGLHS(tmpi, 2:end)); % try again, exit may have been local minimum 
        end
    else     % fitrt == 'fullCDFllh', % always a 2 parameter fit
        initparams=[meanval log(std(XVALS)/3)];  %  guessed matchpoint and sigma
        optima = fminsearch(@fullCDFllh, initparams);
        if max(LOGLHS(:,1)-LOGLHS(end,1) > 0.1)
            [tmp, tmpi] = max(LOGLHS(:,1));
            optima = fminsearch(@fullCDFllh, LOGLHS(tmpi, 2:end)); % try again, exit may have been local minimum 
        end
    end
    bestloglh = LOGLHS(end,1); % best log likelihood and associated parameter values for this condition

%% add extra off-optimal error calculations to help initial estimates of s.e.s; 1D XCs of LOGLHS values via 2D max point:
    % New: base spacing of parameter value xcs on sd of presented
    % stimuli rather than initial estimates as in older versions of psyfit
    offsetmultipliers = [-4 -2 -1 -.5 -.25 0 .25 .5 1 2 4 ]; % variable scales, one should capture cutoff well
    if all(fitrt  == 'fullCDFllh') % use p1xcvals = mean +- multiplier*sigma of presented vals that gave nonuniform responses; OK to be -ve
        p1xcvals = optima(1) + offsetmultipliers.*sigmastimvals{ic};
    else % in POWER_OF_X mode, keep it positive
        p1xcvals = optima(1).*(1 + sigmastimvals{ic}/optima(1)).^ offsetmultipliers ;   % mean * (1 + sigma/mean)^ k
    end
    for k=1:length(p1xcvals)
        errparams=[p1xcvals(k) optima(2:end)]; % when called with 1 parameter value, fitrt adopts DEFAULTEXP
        eval([fitrt,'(errparams);']); % make one XC of 2D parameter space, adding to LOGLHS list of likelihoods and parameter values
    end
    lhvalsp1 = exp(LOGLHS( (end-length(p1xcvals)+1):end,1));
    statstruct{ic}.prelimsigmas(1) = sqrt(var( p1xcvals,lhvalsp1 )); % s.d. of xc of likelihood via mode (at optimum p2)
%{
    NOTES: 
The prelimsigmas, used to form the 2D lattice below, are from the XC of
    likelihood thru the optimum in parameter space. They underestimate
    the sigmas of the marginals if the parameters are strongly correlated.
    By design, the parameters are independent enough so this works.
p1xcvals is itself centrally peaked, not uniformly spaced. SD
    of XC of likelihood underestimates uncertainty for that reason
    (acting like a centrally peaked prior); also liable to underestimate
    by truncation, and (at the opposite extreme) by too coarse sampling.
    Better: estimate curvature of XC of log likelihood from coefts of a
    quadratic fit to the highest triplet in each XC. That avoids sampling
    wobble in the marginals as the lattice points hit or miss the peak
    ridge; and neglecting parameter dependence is a small cost in this
    context, for just the prelimsigmas...But if pdf is heavy-tailed, SD
    might exceed estimate based on quadratic fit to the highest triplets.
    So: take the max (though they mostly seem to agree well):
%}    
    [maxllh, maxllhi] = max(log(lhvalsp1));
    maxllhi = median([2, maxllhi, length(p1xcvals)-1]); % so triplet is definable
    xcllhtripleti = [maxllhi-1:maxllhi+1]; 
    xcllhtriplet = log( lhvalsp1(maxllhi-1:maxllhi+1) )'; % interpolate in margllh from minimal triplet, ignoring tails, for better agreement with s.e.m. of MAP estimates
    quadcoefts = polyfit(p1xcvals(xcllhtripleti), xcllhtriplet ,  2); % all condition-dependent
    statstruct{ic}.sigmaquad(1) = sqrt(-.5./quadcoefts(1)); % approx = sqrt(var(platticevals(2,:), gaussvals(2,:))) etc

    if TWOPARAM=='Y' % two parameter fit 
        p2xcvals=[ -1 -.3 -0.10 -0.04 -0.02 0 +.02 +0.04 0.1 0.3 1] + optima(2); % +, never *, because optima(2) is always a log; scale increment by sigma presented?
        for j=1:length(p2xcvals)
            errparams = [optima(1) p2xcvals(j)];
            eval([fitrt,'(errparams);']);  % make one XC of 2D parameter space, adding to LOGLHS list of likelihoods and parameter values
        end
        lhvalsp2 = exp(LOGLHS( (end-length(p2xcvals)+1):end, 1));
        statstruct{ic}.prelimsigmas(2) = sqrt(var( p2xcvals,lhvalsp2 )); % s.d. of xc of likelihood via mode (at optimum p1)
        [maxllh, maxllhi] = max(log(lhvalsp2));
        maxllhi = median([2, maxllhi, length(p2xcvals)-1]); % so triplet is definable
        xcllhtripleti = [maxllhi-1:maxllhi+1]; 
        xcllhtriplet = log( lhvalsp1(maxllhi-1:maxllhi+1) )'; % interpolate in margllh from minimal triplet, ignoring tails, for better agreement with s.e.m. of MAP estimates
        quadcoefts = polyfit(p2xcvals(xcllhtripleti), xcllhtriplet ,  2); % all condition-dependent
        statstruct{ic}.sigmaquad(2) = sqrt(-.5./quadcoefts(1)); % approx = sqrt(var(platticevals(2,:), gaussvals(2,:))) etc
    end
    statstruct{ic}.prelimsigmas = PRELIMSIGMASCALER*max(statstruct{ic}.prelimsigmas, statstruct{ic}.sigmaquad);
    LOGLHS(:,1)=LOGLHS(:,1)-bestloglh(1); % float vals to zero max, correcting by the max found during iteration
    statstruct{ic}.optima=optima;   % slope parameter is represented as log; this is just to pass the value in statstruct{ic}

%% ************** evaluate llh for 2d lattice of param values, to refine and check 'integrated' (marginal) standard errors
    errvals1= -8:.8 : 8; % was -5:.5:5, gave many warnings...
    % Construct lattice using 'isolated' sigmas from XC of llh; need large multiples of this underestimated sigma
    % sample of 20 is enough to make truncation and sampling errors minimal, shown by
    % agreement with quadratic fit.
    % A non-uniform spacing here would amount to assuming a gently peaked (and
    % post-hoc appropriate) prior in integrating to form the
    % marginals...dodgy, best avoid.
    p1vals=errvals1.*statstruct{ic}.prelimsigmas(1) + optima(1);
    platticevals = p1vals;
    if TWOPARAM~='N' % two parameter fit
        if all(fitrt == 'fullCDFllh')
            errvals2 = -8:.8: 8; % length 21, lengths must match to param 1
        else
            errvals2 = -8:.8: 8; % or sample more finely near center, like a 'post hoc prior'?? Discourages low exponents/shallow slopes, which bias estimate for parameter 1. Dodgy
        end
        p2vals=errvals2.* statstruct{ic}.prelimsigmas(2) + optima(2);  % slope parameter 2 is natural log here, consistent with uniform prior in log
        platticevals = [p1vals; p2vals]; % 2 by nvals
        for k=1:length(errvals1)
            for l=1:length(errvals2)
                errparams=[p1vals(k) p2vals(l)];
                eval([fitrt,'(errparams);']);
                errllh(k,l)=LOGLHS(length(LOGLHS),1);
            end
        end

        errllh=errllh- bestloglh(1); % normalize the likelihood matrix by correcting log by bestloglh, so max llh = 0 and max lh = 1;
        errlh=exp(errllh);  % matrix; linear likelihoods for lattice points, relative to peak of 1
        %{
        Next construct marglh, the marginal absolute probability of the
        data under the model given each sampled value of each parameter
        in turn. This is tricky, requiring averaging over the free
        parameter 2 to obtain the marginal for parameter 1. The average
        is approximated by averaging over the discrete values in the
        sampled lattice. The average should not be a simple unweighted
        one, but should be weighted by the relative probability of the
        different sampled values of the free parameter. For a flat or
        nearly flat prior, that probability is proportional to the
        probability of the data , as a fraction of the sum of the data
        probability for all values of the free parameter when the fixed
        parameter is held at x ('empirical Bayes'). The 'empirical Bayes'
        estimate of the probability of the data at param1 = x is
        the mean for all y at fixed x of
            P(x,y) = p(data | param1 = x, param2 = y )
          with weights given by P(x,y) divided by the sum of
        P(x,y) for all y.
        This 'empirical Bayes' weighting tightens the marginal likelihood
        distributions and makes s.e.m.s agree well with empirical
        reproducibility of the MAP estimates. More extreme weightings
        than this give about the same results.
        %}
        %        errlh = errlh.*exp(bestloglh); % errlh is p(data given each pair of parameter values in lattice), or P(x,y) above
        statstruct{ic}.marglh = [sum(errlh'.^2)./(sum(errlh') + realmin); sum(errlh.^2)./(sum(errlh) + realmin)]; %  [ p(data | x ); p(data | y )] as 2 by nvals;
        %        realmin is around 10^-308, system dependent, to avoid /0
        %        which seldom arises
        if MAKEPLOTS>1 % show lh as f(param 1, param2)
            figure(ic);
            contour(p1vals,exp(p2vals),errllh',[-.5 -1 -2 -4 -10 -40]),title('contours of (errllh) at ln lh = [-.5 -1 -2 -4 -10 -40]');
            grid;
            shg;
            %            disp('LLH contour figure created; dock the figure window if you want to keep figures visible; press any key to continue');
            %            pause;
            %            mesh(p1vals,exp(p2vals),errlh'),title('mesh(errlh)')
            %            disp('LLH surface figure created; press any key to continue');
            %            pause;
        end
    else %  single parameter fit
        l = 1;
        for k=1:length(errvals1)
            errparams = p1vals(k);
            eval([fitrt,'(errparams);']);
            errllh(k,1)=LOGLHS(length(LOGLHS),1)- bestloglh(1); % normalize the likelihood matrix by correcting log
        end
		errlh = exp(errllh);
        statstruct{ic}.marglh = errlh'; % row vector of data probabilities,  for quadratic fit below
    end % of 2nd "if TWOPARAM... else..  block, for creating llh for lattice points
    LOGLHS((end-k*l+1):end,1) = LOGLHS((end-k*l+1):end,1) - bestloglh(1); % normalize LOGLHS values to about zero max,

%% ***************get mean and sd of marginal pdfs, and fit 1D quadratics to the loge marginal  pdfs:
    np = size(platticevals,1);
    meanquad = NaN(1,np); % will be mean from quadratic fit to ln likelihood (for lattice values only)
    for ip = 1: np
        lh = statstruct{ic}.marglh(ip,:); % probability of the particular response sequence, for given value of parameter ip, other parameter unconstrained.
        % fit the peak triplet with a quadratic for interpolation:
        [maxllh, maxllhi] = max(lh);
        maxllhi = median([2, maxllhi, length(lh)-1]); % so triplet is definable
        margllhtripleti = [maxllhi-1:maxllhi+1]; 
        margllhtriplet = log( lh(maxllhi-1:maxllhi+1) ); % interpolate in margllh from minimal triplet, ignoring tails, for better agreement with s.e.m. of MAP estimates
        quadcoefts(ip,:) = polyfit(platticevals(ip,margllhtripleti), margllhtriplet ,  2); % all condition-dependent
        % curvature (1st val) might not be negative at edge: kluge if so
        if quadcoefts(ip,1)>=0
            meanquad(ip) = platticevals(1,ip)*(maxllhi==2) + platticevals(end,ip)*(maxllhi==length(platticevals)-1) % put at edge
            quadcoefts(ip,1) = -eps; % model llh as exponential from edge value
        else
            meanquad(ip) = - 0.5*quadcoefts(ip,2)./ quadcoefts(ip,1); % "-b/2a" % mean is least squares estimate of best paramval for any distribution
        end
        % Natural log likelihood drops by .5 at 1 sd, this is quadcoefts(ip,1)*(pvals-mu).^2, so at +-sigma, (pvals-mu).^2 = -.5/quadcoefts(ip,1)
        quadvals(ip,:) = polyval(quadcoefts(ip,:),platticevals(ip,:)); % the quadratic fit to log(marglh) around the peak
        gaussvals(ip,:) = exp(quadvals(ip,:)); % Gaussian fit to statstruct{ic}.marglh, based on curvature near peak 

        % quantiles by interpolation...cumtrapz for integral of lh; ensure
        % distinctness of values in lhint for interp1, as in fastEstimate()
        quantileLevels = [.005 .025 .05 .5 .95 .975 .995];

        lh = statstruct{ic}.marglh(ip,:)./sum(statstruct{ic}.marglh(ip,:)); % normalization just to get nonzero reports in warning. Use prior cubic interp in log(marglh)?
        if ((lh(1) > max(lh)*10^-2) || (lh(end) > max(lh)*10^-2)) % ^-4?? dm ^-2.35 is roughly 99.9% confidence interval, no? should be sufficient -ev
            fprintf('\nWarning: edges not much less likely than max (param %d): \n\tmaxlh = %f, edges: [%f, %f]\n\t If extreme quantiles are important,\n\ttry a more extended lattice range.\n\t', ip, max(lh), lh(1), lh(end));
            %            keyboard;
        end  % Or try redo from line 425 or so with spacing based on presented values: new function, or new mode...
        %        Call recursively in alternate mode, flagging that to avoid further tries??
        interpx = linspace(platticevals(ip,1), platticevals(ip,end), 10*(size(platticevals,2)-1)+1); 
        interplh = exp(  interp1(platticevals(ip,:),log(lh), interpx, 'spline')  ); % use cubic in log for interp
        lhint = cumtrapz(interpx,interplh);   % rel. fraction of area under lh distn <= x, by approx integration
        %         using trapezoid rule
        lhintN = lhint./max(lhint); % fraction of area under lh distn <= x, by approx integration
        [cumlk, tmpi] = unique(lhintN);
        statstruct{ic}.quantiles(ip,:) =  interp1(cumlk,interpx(tmpi),quantileLevels,'spline');
    end
    % Similar result with [p,s,mu] = polyfit(platticevals(ip,okmargidx), margllhok,  2); quadvals = polyval(p,(platticevals(ip,:)-mu(1))./mu(2)   )
    statstruct{ic}.meanquad = meanquad;  % ic x nlatticevals x nparams
    statstruct{ic}.sigmaquad = sqrt(-.5./quadcoefts(:,1))'; % approx = sqrt(var(platticevals(2,:), gaussvals(2,:))) etc
    statstruct{ic}.meanlatticemarg = sum( statstruct{ic}.marglh.* platticevals,2 )'./(sum(statstruct{ic}.marglh,2)');  % row, 1 or 2
    ms_mlh = sum( statstruct{ic}.marglh.* (platticevals.^2),2)./sum( statstruct{ic}.marglh,2); % row, size nparams
    statstruct{ic}.sigmalatticemarg = sqrt(ms_mlh - statstruct{ic}.meanlatticemarg( :).^2)'; % row vector

%% convert 2nd parameter and s.e. from log value, put in statstruct{ic}
    statstruct{ic}.platticevals(1,:) = platticevals(1,:);
    if np > 1 % we have a second free parameter
        statstruct{ic}.optima(2) = exp(statstruct{ic}.optima(2)); %   convert back fminsearch 2D maxlh estimate from the log form passed to fitrt
        statstruct{ic}.platticevals(2,:) = exp(platticevals(2,:));
        statstruct{ic}.meanquad(2) = exp(meanquad(2));
        statstruct{ic}.sigmaquad(  2) = 0.5*(exp(meanquad(2)+ sqrt(-.5./quadcoefts(2,1))) - exp(meanquad(2)- sqrt(-.5./quadcoefts(2,1)))); % piecewise linear approx to log
        statstruct{ic}.sigmalatticemarg(  2) = 0.5*(exp(statstruct{ic}.meanlatticemarg(2)+ statstruct{ic}.sigmalatticemarg(2))  - exp(statstruct{ic}.meanlatticemarg(2)- statstruct{ic}.sigmalatticemarg(2)  ));
        statstruct{ic}.meanlatticemarg(2) =  exp(statstruct{ic}.meanlatticemarg(2));
        statstruct{ic}.quantiles(2,:) =  exp(statstruct{ic}.quantiles(2,:));
    end

%% ******************** derive outputs, make plots if MAKEPLOTS>1
    % Derive fits (actual and predicted proportions correct):
    USEQUAD = 1; % rather than raw sd of lattice marginal pdf values
    USEMAP = 1; % report 2D max a posteriori estimates as bestparams...
% ...this choice adopts  the MAP estimates as ‘best’ and uses them for the
% fitted curve, but still uses the marginals for the standard errors and
% quantiles of the parameter values. With the ‘marginals’ option, not only
% is the fit poorer but the ‘best’ threshold and ‘best’ sigma are not
% generally *jointly* best. But the marginal estimates are better for
% minimizing the cost of errors in the estimates where cost is graded with
% error size. Usually the data constrain the estimates well enough that it
% doesn’t make much difference.
    if USEQUAD % alternative is raw sd of lattice marginal pdf values, legit if lattice is uniformly spaced
        bestparams(ic, :) =[statstruct{ic}.meanquad  statstruct{ic}.sigmaquad ];
        statstruct{ic}.bestparams = statstruct{ic}.meanquad;
        %  ...use original optima (MAP estimates for mean, closer to actual
        %  sample data...but heavy tailed pdfs could favor USEQUAD=0)
    else % don't use quadratic approx to marginals:
        bestparams(ic, :) =[statstruct{ic}.meanlatticemarg  statstruct{ic}.sigmalatticemarg ]; % 4 long vectors for each condition as summary result: mean 1, mean 2, se1, se2; for each condition
        statstruct{ic}.bestparams = statstruct{ic}.meanlatticemarg ;
    end
    if USEMAP % prefer MAP (2D optimum) for best parameter value, keeping sigma estimate
        bestparams(ic, 1:length(statstruct{ic}.optima)) = statstruct{ic}.optima; % preserve sems as  elements 3:4 if present
        statstruct{ic}.bestparams = statstruct{ic}.optima;  % nparams (one or two) elements
    end
    statstruct{ic}.curvex = linspace(xpresented{ic}(1), xpresented{ic}(end), 40); 
    if (length(statstruct{ic}.bestparams) >  1) 
        [tmp ppc{ic}] = eval([fitrt,'([statstruct{ic}.bestparams(1) log(statstruct{ic}.bestparams(2))]);']); % define predictedpc; fitrt uses log transformed value for shape parameter
    else
        [tmp ppc{ic}] = eval([fitrt,'(statstruct{ic}.bestparams(1))']); 
    end
    statstruct{ic}.fits = [xpresented{ic} Npresented{ic} opc{ic} ppc{ic}]; % for return
    statstruct{ic}.probdata = sum(sum(errlh.^2))./sum(sum(errlh)).*exp(bestloglh);  % prob of all data under this model, absolute
    statstruct{ic}.marglh = statstruct{ic}.marglh.*exp(bestloglh);  % make marglh also an absolute average probability
    lapseprob= 0.025; % should be global for consistency with fitrt 
    if fitrt=='halfCDFllh'
        if TWOPARAM=='Y'
            n = statstruct{ic}.bestparams(2) ; %exponent on x; only if forced choice
        else
            n=DEFAULTEXP;
        end        
        xn=(1/sqrt(2.0))./(statstruct{ic}.bestparams(1)./statstruct{ic}.curvex).^n;
        statstruct{ic}.curvey = 1/NCHOICE + (1-lapseprob)*(1-1/NCHOICE)*erf(xn);
    else % assume usual fullCDFllh
        xdev=(statstruct{ic}.curvex-bestparams(ic,1))./(sqrt(2.0).* statstruct{ic}.bestparams(2));
        if NCHOICE>1 % nafc detection, treat lapses asymmetrically
            statstruct{ic}.curvey = 1/NCHOICE +(1-1/NCHOICE)*(0.5.*(1-lapseprob).*(1+erf(xdev))); % stay on 50% for chance,
        elseif NCHOICE==0  % comparative judgment
            fingerrprob= 0.02; % to allow for lapses with extreme stimuli; check same as in fitrt
            statstruct{ic}.curvey=0.5.*(1+(1-fingerrprob).*erf(xdev)); % ppc is a normal integral: see halfCDFllh.m. Treatment of fingerprob symmetrical, could be improved
        else  % NCHOICE = 1, yes/no detection 
            statstruct{ic}.curvey = 0.5.*(1-lapseprob).*(1+erf(xdev)); % stay on 50% for chance, 
        end
    end
    LOGLHS=LOGLHS(1:length(LOGLHS)-1,:); % for plotting, trim LOGLHS from this last fitrt call, because unnormalized; also trim statstruct:
    statstruct{ic} = rmfield(statstruct{ic}, {'prelimsigmas'}); % wanted only for intermediate computations
    if MAKEPLOTS>1 % visually check fit of quadratic to log likelihoods
        % Main deviation is that predLOGLHSorth assuming orthogonality (via
        % precs_orth) overestimates the probability of the data for the XC
        % scans thru the 2D max point; this is expected
        figure(nconds+ic)
        for ip = 1: np
            subplot(np,2,2*ip-1);
            plot(statstruct{ic}.platticevals(ip,:), statstruct{ic}.marglh(ip,:),'--*') % complete marglh
            hold on
            grid
            plot(statstruct{ic}.platticevals(ip,:),gaussvals(ip,:),'-o')
            if ip ==2
                paramstring = 'slope or exponent';
            else
                paramstring = 'parameter 1';
            end
            title(['fit of gaussian (o) to marginal LH (*) vs. ', paramstring ]);
            subplot(np,2,2*ip)
            plot(statstruct{ic}.platticevals(ip,:), log(statstruct{ic}.marglh(ip,:)),'--*');
            hold on
            grid
            plot(statstruct{ic}.platticevals(ip,:),quadvals(ip,:),'-o');
            title(['fit of parabola (o) to ln marginal LH (*) vs. ', paramstring ] );
        end
        shg
        hold off
    end
end  % of ic=1:nconds loop to define bestparams matrix ( 1 row per condition), plus statstruct structure and fits for each condition


%% ******* clean up, announce best fit values, plot fits for all conditions
if MAKEPLOTS % do this plot by default (MAKEPLOTS nonzero)
    figure
    sym = {'+', 'o', 'x','*'};
    plotsym = {'ro--' 'go--' 'bo--' 'ko--' 'rs--' 'gs--' 'bs--'  'ks--' 'rd--' 'gd--' 'bd--' 'kd--'  };
    for ic = 1:nconds
        legs{ic} = ['Condition ' num2str(ic)];
        h(ic) = plot(statstruct{ic}.curvex, statstruct{ic}.curvey,...
            plotsym{1+mod(ic-1,length(plotsym))}([1 3]), 'LineWidth',2);% continuous curve,no symbols
        hold on;
        for idx = 1:length(xpresented{ic})
            plot(xpresented{ic}(idx),opc{ic}(idx),plotsym{1+mod(ic-1,length(plotsym))}([1 2]), 'MarkerSize', Npresented{ic}(idx)+2  ); % could equally use .fits here
        end
        plot(statstruct{ic}.fits(:,1), statstruct{ic}.fits(:,3),... % opc
        plotsym{1+mod(ic-1,length(plotsym))} ,'MarkerEdgeColor',plotsym{1+mod(ic-1,length(plotsym))}(1), 'MarkerFaceColor', plotsym{1+mod(ic-1,length(plotsym))}(1),'MarkerSize',8)% observed and predicted % correct for each x presented.
        plot(xpresented{ic},opc{ic}',plotsym{1+mod(ic-1,length(plotsym))}([1 3 4]));
    end
    legend(h,legs,'Location', 'SouthEast');
    xlabel('staircased stimulus value');
    ylabel('proportion positive responses');
    shg
    hold off; % to avoid overwriting this plot

    disp('best parameter 1 (and 2 if 2 parameters varied), followed by sem1 (and sem2 if applicable), in 2 0r 4 cols, 1 row per condition:');
    disp(bestparams);
end

if length(statstruct) ==1
    statstruct = statstruct{1};
end



%%
function [llh, ppc] = halfCDFllh(paramvector)
%Function llh = llh(paramvector) returns minus the natural log of the
% likelihood of frequency-of-seeing data from globals nc, ntt
% for given sigma and loge(exponent) on global xvals; called by psyfit.m
% Sigma=paramvector(1),loge(exponent)=paramvector(2).
% Uses reciprocal of xval if data suggest that % correct decreases with increasing xval. This rather than a sign change
% on xval seems right for the fc case: there has to be a true zero, on x or else on its reciprocal.
% Also appends llh and assocd param vals to LOGLHS list
% LOGLHS ends up with dimensions ITERATIONs x 3 (make global to access).
global TWOPARAM DEFAULTEXP NCHOICE;
global XVALS NTT NC LOGLHS ITERATION;
PRIORMAXLOGEXP = log(2.5); PRIORSIGMA = 0.4; % gentle Gaussian default prior on log(exponent), lh = x 0.8 if exponent is 1.4 x max or .7 x max; could request user input for prior...
sigma = max(paramvector(1), eps);
if TWOPARAM=='Y'
    n=exp(paramvector(2)); %exponent on x; only if forced choice
else
    n=DEFAULTEXP; % fixed global exponent for single parameter fit, varying threshold only
end
paramvector(2)=log(n); % will appear in column 3 of LOGLHS
logprior = -(log(n)-PRIORMAXLOGEXP)^2/(2*PRIORSIGMA^2); % flattish Gaussian in exponent, flat for threshold...
if ITERATION == 1
    meanval = (XVALS'*(NTT))./sum(NTT);
    meancval = (XVALS'*(NC))./sum(NC);
    if meancval <= meanval   % correctness-weighted mean minus unweighted mean; use reciprocal if negative
        disp('data are not clearly monotonic increasing; consider taking reciprocals of stimulus values');
    end
    LOGLHS=[];
end
% now keep the 84% point as a pivot while varying n
xn=(1/sqrt(2.0)).*(XVALS./sigma).^n;
% note: matlab erf is not a normal integral, despite what they say; next line
% makes ppc = (lapseprob * 0.84) if xn = .707, or if xval = sigma above. So ppc
% is a normal integral*(1-lapseprob), of (xval/sigma)^n, with usual sigma in xval.
lapseprob= 0.025; % to allow for lapses, including finger errors with extreme stimuli, treated as asymmetrical
% Fitted p ranges from chance to just over 1-lapseprob 
ppc = 1/NCHOICE + (1-lapseprob)*(1-1/NCHOICE)*erf(xn);

llh=-logprior-sum(NC.*log(ppc+eps)+(NTT-NC).*log(1-ppc+eps)); % actually = -log(llh) ; why? for minimization, historic?
LOGLHS=[LOGLHS; -llh paramvector(1) paramvector(2)]; % append real loglh, etc
ITERATION = ITERATION+1;

%%
function [llh, ppc] = fullCDFllh(paramvector)
%function llh = fullCDFllh(paramvector) returns minus the natural log of the
% likelihood of frequency-of-seeing data in globals nc, ntt
% for given 50% and loge(sigma); called by psyfit.m
% 50% point (subjective match,etc.)=paramvector(1); for yes/no, not f/choice.
% Also appends llh and assocd param vals to loglhs list
% loglhs ends up with dimensions iterations x 3 (make global to access).
% Second parameter (psychometric width) is passed as a natural log, and
% exponentiated here, to avoid going negative during iteration.
global NCHOICE;
global XVALS NTT NC LOGLHS ITERATION  ;
PRIORMAXLOGSIGMA = 0.5*log(var(XVALS, NTT)/2); PRIORSIGMA = 1; USEPRIOR = 0; % gentle Gaussian default prior on natural log of psychometric sigma. BUGFIX: was std(XVALS.*NTT)...
% most likely sigma is set to half the sigma of presented stimuli, uncertainty around that is about a factor of two for a reduction of prior to 80%
% could request user input for prior...
matchval=paramvector(1);
sigma=exp(paramvector(2)); % sigma for psychometric function
if USEPRIOR
    logprior = -(log(sigma)-PRIORMAXLOGSIGMA)^2/(2*PRIORSIGMA^2); % flattish Gaussian  prior for psychometric sigma, no (or flat) prior for threshold...
else 
    logprior = 0;
end
if ITERATION == 1
    LOGLHS=[];
    meanval = (XVALS'*(NTT))./sum(NTT);
    meancval = (XVALS'*(NC))./sum(NC);
    if meancval <= meanval   % correctness-weighted mean minus unweighted mean
        warn('data are not clearly monotonic increasing; consider sign swap for stimulus values');
    end
end
fingerrprob = 0.02; % .97; to allow for finger errors with extreme stimuli in comparative judgment
lapseprob= 0.025; % to allow for detection lapses with extreme stimuli
xdev= (XVALS-matchval)./(sqrt(2.0).*sigma);
if NCHOICE>1 % nafc detection, treat lapses asymmetrically
    ppc = 1/NCHOICE +(1-1/NCHOICE)*(0.5.*(1-lapseprob).*(1+erf(xdev))); % stay on 50% for chance, 
elseif NCHOICE==0  % comparative judgment
    ppc=0.5.*(1+(1-fingerrprob).*erf(xdev)); % ppc is a normal integral: see halfCDFllh.m. Treatment of fingerprob symmetrical, could be improved
else  % NCHOICE = 1, yes/no detection 
    ppc = 0.5.*(1-lapseprob).*(1+erf(xdev)); % stay on 50% for chance, 
end
llh=-logprior-sum(NC.*log(ppc+eps)+(NTT-NC).*log(1-ppc+eps)); % actually = -llh
LOGLHS(ITERATION,1)=-llh;           % real loglh
LOGLHS(ITERATION,2)=paramvector(1);
LOGLHS(ITERATION,3)=paramvector(2);
ITERATION = ITERATION+1;


% DETAILS:
% It is assumed that response probability, nc/ntt, increases with xval.
% S.e.s of parameters are from binomial variability only, so 
% session-session variability will be greater, but these s.e.s are
% appropriate for comparing conditions randomly interleaved with sessions.
%
% For other than the 2afc POWER_OF_X OPTION, psyfit fits a full ogive,
% with generally non-zero intercept; in nAFC cases, a compressed ogive is
% fit going to chance at x = -Inf; see above. For the POWER_OF_X option,
% halfCDFllh is called with an exponent on x as the second parameter. Otherwise
% fullCDFllh is called, and the second parameter is the s.d. of the Gaussian
% corresponding to the fitted cumulative normal.
%
% The second parameter (exponent for halfCDFllh, curve slope for
% fullCDFllh) is passed as its logarithm, partly to prevent sign reversals.
% The exponent is poorly constrained in runs of 100 or so, and the best fit
% can be very high, and averages of it can be biased.
% Standard errors are found by fitting llh with quadratic in first
% parameter and in log-transformed 2nd parameter. Interpretation is a bit
% uncertain in both cases once the s.e. becomes comparable with the optimum
% value..NaN is returned as the s.e. if no useful estimate can be made; in
% that case try repeating the analysis with TWOPARAM='N'.
%  Quadratic description of log likelihood is illustrated if
%  MAKEPLOTS > 1
% It often fails in the tails, but s.e. estimates don't depend on tails
% much. 
% The variation of likelihood over the 2D parameter surface is assessed.
% S.e.s are then from the marginal distribution of likelihood; this allows
% for some redundancy between parameters, since the uncertainty in each
% parameter is evaluated on the assumption of complete uncertainty about
% the other. If ORTHOG= 'Y', independence of the parameters is assumed. For
% speed, it usually is not disastrous to evalue s.e.s 1-dimensionally,
% which is done by evaluating likelihood as f(parameter of interest) with
% the other parameter at its best-fitting value. This is done for each
% parameter in turn after iterations converge, and the results can be
% visualized if MAKEPLOTS > 1. Generally, 1D s.e.s are naturally
% smallest, followed by ORTHOG s.e.s, then s.e.s based on fits allowing
% interaction between parameters, with likelihood matrix-based s.e.s
% largest. The latter are preferable, but slower, and may be inflated by
% allowing interaction which is generally in fact small, but which can be
% accidentally enough to inflate the obtained s.e. Matrix s.e.s with the
% orthogonality constraint are close to quadratic unconstrained ones and to
% simulation-based ones. (Matrix/quadratic s.e. ratio is about 1.2:1, and
% unconstrained:orthogonally constrained ratio is about the same, so these
% non-optimalities cancel). Adoption of the orthogonality constraint would
% usually be OK but can underestimate s.e.s in cases where only the tail or
% part of the frequency of seeing curve is available. 
%{
% Note on BEST STIMULUS PLACEMENT: in fitting a full ogive to y/n data the
% steepest gradient is at the inflection or 50% point, and this is also
% the stimulus value at which the precision in the estimated log
% sensitivity is greatest, since binomial variation in p is also least at
% p - .5. For 2AFC, gradient is least at p = .75, but the shrinkage in 
% binomial variability at increasing p shifts the best point to .83, where
% the product of binomial sd in p with gradient of p with log sensitivity
% is least. In power of x mode, the +1 sigma point would be best if 
% binomial variation were constant, and the optimum p is slightly
% upshifted from 0.83 as discussed by D.M. Green, JASA 1960s.
z = -2:.1:+2
p = (1 +(1 +erf(.707*z))/2)/2  % p vs z 
g = exp(-z.^2/2)  % slope of p vs z, reciprocally affects precision of est.
sdlogsens = sqrt(p.*(1-p))./g % sd of p, divided by slope for sd of z
[minsd, bestz] = find( sdlogsens ==min(sdlogsens)) % sweet spot 
p(bestz)
This sort of analysis yields about p = .7 for 4AFC,  and p = 0.85 or so for
2AFC. This assumes a model for the psychometric function where proportion
correct is a guessing-corrected cumulative normal in *log* stimulus
value (because only then does the 'sweet spot' minimize error in the 
log of the threshold estimate, rather than in the threshold value itself).

For 2AFC, POWER_OF_X mode fits p with a cumulative normal function of
linear stimulus value. For that mode, my analysis agrees with Laming and
Marsh (and some earlier papers) in suggesting a sweet spot at a much higher
optimal p, about 0.94. Nobody seems to actually track these high
probabilities, but I believe they really are more efficient. I like 0.85 or
so for 2AFC well enough, in rough conformity with the log model, even though
the linear model ideally calls for a more extreme target probability.
DRJ Laming and D. March, Perception and Psychophysics 1988...
%} 
%
%
% RECOMMENDATIONS:
%? USEQUAD = 1; else direct marginal pdf sd and mean (latticemarg fields of
% statstruct) areinfluenced by truncation and sampling bias. ??
%(MAKEPLOTS = 1; % set to > 1 to obtain intermediate plots, 0 for no plot
%(TWOPARAM='Y'; % set to 'Y' to obtain 2 parameter f/c fits;'N' for one parameter
%(TWOPARAM can be overwritten by user dialog in MAKEPLOTS (default) mode.
%(TWOPARAM and MAKEPLOTS can be overwritten by supplied arguments...
%(DEFAULTEXP=2.5; % used only if varyexp='N' in POWER_OF_X mode
%(ORTHOG='N'; % don't assume orthogonality: in a few cases it fails
% The quadratic sigmas may have the wrong sign in an unconstrained fit;
% the script then resorts to an orthogonal-constrained fit. (This can happen
% when there are local minima, or as an accidental result of placement of
% iteration values, creating a large apparent tradeoff between param vals.)
% If this should ever be insufficient to deliver an s.e. estimate, a
% possible recourse is to run the offending condition with a fixed exponent.

%{
*****************CHANGES:

psyfit Aug 2012, differs from prior version in that marglh element of
statstruct is the probability of the data for the sampled values of each
model parameter, using empirical Bayes weighting over the free parameter.
POWER_OF_X mode is usable in any FC experiment (not just 2AFC as before).
Thus with nchoice == 4, and a non-zero supplied value for exponent,
psyfit uses ‘POWER_OF_X’ mode to fit the z >= 0 half of a normal CDF
to the probability of seeing for linear stimulus values which MUST all be positive.
But as always when only half the CDF is fit, the generally needed
inflection in the theoretical curve is lacking unless the  stimulus value
is taken to an exponent > 1. In psyfit.m, this was possible only for 2AFC
since there was a rationale for it (Gaussian noise added to a power law
function of stimulus value: Kontsevich and Tyler). Psyfit now does
the same for other nchoice values, eg nchoice = 4, simply by scaling the
theoretical curve for chance at stimulus value zero.

CHANGES, 2014: Peak of marginal likelihoods uses only the minimum 3 values.
Cubic spline interpolation prior to cumtrapz for more accurate integration.
When marginal peaks at edge, positive curvature of llh distribution is 
ignored, modeling llh as linear (lh as exponentially decreasing from edge).
Sept 2014: use cross-sections thru mode to determine lattice spacing, from
peak curvature or sd, whichever is greater.
Use marginal distributions rather than quadfit2D (discarded) to find sems
and quantiles for parameter values. 

% ***********POSSSIBLE ENHANCEMENTS:
% 1. 75% (halfway) thresholds rather than 84% (1 x sigma)? Just x .6745...
% 2. replace globals with passed parameters (on aesthetic grounds only).
% 3. allow other cdfs than Gaussian (e.g. from FAST?)...not critical.
% 4. line 553 Call recursively in alternate mode, flagging that to avoid further tries??
% 5. line 425 switch for POWEROFX for add vs scale
% 6. line 438 p2xcvals=[ -1 -.3 -0.10 -0.02 0 +.02 +0.04 0.1 0.3 1] + optima(2); scale increment by sigma presented?
% 7. "Prior' on sigma could be introduced to aid convergence, and if
% necessary removed from final marginal pdf as f(sigma). Fitrt needs access
% to parameters for prior in that case (globals or arguments).
% 8. Incorporate reflected/inverted stimulus values into output fits where
% data have p decreasing...reflection for fullCDF, inverted in POWER_OF_X.
% 9. Avoid spline interp on theoretical curve, evaluate at curvex values directly
% 10. Reflections and inversions, assign matchval = WHICHSIGN*paramvector
% XVALS reflected, restore for plotting.
%  


%% Worries and reassurances:
% Stimulus values expressed as logarithms should only be used in full ogive
% mode, but this mode can also be appropriate when linear stimulus values
% all fall far from zero, as in many matching experiments. The full ogive
% is constrained to chance at zero stimulus value. Generally any ogive will
% do for a useful description and parameterization, but detection data
% often deviate from odd symmetry around the middle of the ogive, in a way
% that the POWER_OF_X mode can pick up. In full ogive mode, the normal cff
% of log stimulus magnitude has the opposite asymmetry, so there is some
% small incentive to consider alternatives like a Weibull shape that fit
% typical detection data a little better. Not a big deal in my view.
%
% Notches occur in predloglhs (but not in predloglhsorth) when threshold
% parameter is off-optimum, and slope parameter is scanned through optimum:
% These are expected when the precs coefficients are inappropriate in
% that part of parameter space. Where the interaction term gets too much
% weight, a symmetrical notch is the necessary result.
%
%
%}