function [ThreshDiffLkhd, WidthDiffLkhd, BayesFactor2params, bestparams, statstructs, bestparamsmerged, statstructmerged ] = threshdifflkhd(twoconditiondata, nchoice)
% function [ThreshDiffLkhd, WidthDiffLkhd, BayesFactor2params, bestparams, statstructs, bestparamsmerged, statstructmerged ] = threshdifflkhd(twoconditiondata, nchoice)
% Run psyfit on data from two conditions, get Bayes Factors for models
% constrained to equality of psychometric function parameter 1 (translation)
% or parameter 2 (width or slope), or both.
% Uses psyfit, with full normal cdf psychometric function and with scaling
% for 'guessing' if required by nchoice. 
% Set nchoice to 1 for yes/no or pse experiments.
% This treatment is only loosely Bayesian, with implicitly flat prior.
%
% Results support simpler computations using parameter s.e.m.s.  
% Trials suggest these procedures give essentially the same answer as using 
% the standard error of each threshold from psyfit...the pdf of delta has the 
% same standard deviation as the quadratic sum of the s.e.m.s of the
% thresholds estimated by psyfit for each condition. So this is mainly
% reassuring rather than a substantial improvement in method.

% Form 2 separate statstructs for independent fits to the two conditions:
[bestparams, statstructs ] = ...
psyfit(twoconditiondata, nchoice, 0, 0 );

% Fit merged data with a single pair of parameter values
mergeddata = twoconditiondata; 
% now merge data by resetting condition indices to 1:
mergeddata((mergeddata(:,4) ==2),4) = 1;
[bestparamsmerged, statstructmerged ] = ...
psyfit(mergeddata, nchoice, 0, 0 );

% 1. Test for a difference in either parameter (or both):
% The probability of the data under a model, assessed with a
%likelihood-weighted average over model parameter values, is
%available in statstruct.probdata. For the unconstrained model, where
%threshold values can differ for the two conditions, the probability of the
%combined data, as a function of both the threshold estimates t1 and t2 for
%the two conditions, can be obtained by multiplication, as
%statstruct1.probdata*statstruct2.probdata from separate analyses of data
%sets D2 and D1. For the constrained model, psyfit above applied to 
%the merged data yields statstructmerged.probdata. The ratio of this
%to the product is the relative probability of the data under the
%constraint vs without the constraint. This is probably not much different
%from the corresponding ratio of probabilities based on
%the MAP parameter values for each of the 3 cases, without integrating
% over the likely range of parameter values; ratios of peaks may be fairly
%close to ratios of volumes. But this is better and done already in psyfit.

BayesFactor2params = statstructmerged.probdata./(statstructs{1}.probdata*statstructs{2}.probdata);

% THis is the relative likelihood of the data assuming common values 
% for the psychometric function parameters (threshold and width), as 
% compared with the assumption of independent optimal values for each
% condition.

% Test for a difference in inflection point of ogive (pse) with slope equated:
% This is surprisingly hard to do correctly.
% The main difficulties dealt with here are
% (1) to integrate appropriately over the common but uncertain slope
% or width parameter, and 
% (2) to integrate over uncertainty in the threshold parameter itself, 
% in assessing the likelihood of a given difference in threshold.
% See comparingthreshold.docx for more detailed explanation...
%
% Note: A simpler test is to consider the probability of each condition's
% data and of the merged data assuming the MAP parameter values in each 
% case. The ratio of the merged probability to the product of the component
% probabilities is analogous to the ratio computed below, and often close. 
% Here we do take the average over uncertain parameter values.

t1array = statstructs{1}.platticevals(1,:);  % row vector of t for cond. 1
t2array = statstructs{2}.platticevals(1,:);  % row vector for cond. 2
ranget1vals = t1array(end) - t1array(1); 
newt1vals = [statstructs{1}.platticevals(1,:)  statstructs{1}.platticevals(1,1:end-1)+ 0.5*diff(statstructs{1}.platticevals(1,:))]; % add intermediate vals 
newt1vals = (sort(newt1vals)); % sort in order, row vector for interpolation
delta = linspace(-ranget1vals, ranget1vals, 2*length(newt1vals))';  % column vector

% Get prob2 for newt1vals + various deltas, and
% integrate prob2(t1 +delta)*prob1(t1) with respect to t1...
% YI = INTERP1(X,Y,XI)
% Linear platticevals insufficient for interpolation, hence negative
% interpolated probabilities. Interpolation in log works OK

prob1 = statstructs{1}.marglh(1,:);
prob2 = statstructs{2}.marglh(1,:);
for idx = 1:size(delta)
    log10prob1new(:,idx) = interp1(t1array', log10(max(prob2',realmin))', newt1vals,'spline'); % Condition1, prob1new is Nnewt1vals x Ndeltas, row m is constant at prob1(newt1vals(m)), no variation with delta
    log10prob2new(:,idx) = interp1(t2array', log10(max(prob2',realmin)), newt1vals'+delta(idx),'spline'); % Condition2, prob2new(m,n) is prob2 interpd at (newt1vals(m) + delta(n))
end
prob1new = 10.^log10prob1new;
prob2new = 10.^log10prob2new;
prob1and2new =  prob1new.*prob2new; % pointwise, not matrix product...prob of all data together for t1 = newt1vals(m), t2 = t1 + delta(n)

ThreshDiffLkhd= [delta mean(prob1and2new,1)'./max(mean(prob1and2new,1))];  
% Row 2: Likelihoods for TRANSLATION of psychometric function between 
% the compared conditions (normalized to peak of 1 in the array)
% Row 1: corresponding 
% translations (delta = difference in the threshold parameter).
% The value for zero delta is a likelihood ratio for the 'same threshold' 
% hypothesis compared with the post-hoc most likely threshold
% difference.
% Could get pdf for delta by normalizing area, and fiducial limits by 
% cumtrapz, but this hardly seems worthwhile since simpler substitutes
% work well.

% Test for a difference in SLOPE of ogive, with inflection point pinned:
% as above but swapping the roles of params 1 and 2 for each case...
% spline here really wants log of slope parameter, for regularity of
% spacing, but this should not be critical...
w1array = statstructs{1}.platticevals(2,:);  % row vector
w2array = statstructs{2}.platticevals(2,:);  % row vector
rangew1vals = w1array(end) - w1array(1); 
neww1vals = [statstructs{1}.platticevals(2,:)  statstructs{1}.platticevals(2,1:end-1)+ 0.5*diff(statstructs{1}.platticevals(2,:))]; % add intermediate vals 
neww1vals = (sort(neww1vals)); % sort in order, row vector
delta = linspace(-rangew1vals, rangew1vals, 2*length(neww1vals))';  % column vector

prob1 = statstructs{1}.marglh(2,:);
prob2 = statstructs{2}.marglh(2,:);
for idx = 1:size(delta)
    log10prob1new(:,idx) = interp1(w1array', log10(max(prob2',realmin))', neww1vals,'spline'); % Condition1, prob1new is Nneww1vals x Ndeltas, row m is constant at prob1(neww1vals(m)), no variation with delta
    log10prob2new(:,idx) = interp1(w2array', log10(max(prob2',realmin)), neww1vals'+delta(idx),'spline'); % Condition2, prob2new(m,n) is prob2 interpd at (neww1vals(m) + delta(n))
end
prob1new = 10.^log10prob1new;
prob2new = 10.^log10prob2new;
prob1and2new =  prob1new.*prob2new; % pointwise, not matrix product...prob of all data together for w1 = neww1vals(m), w2 = w1 + delta(n)

WidthDiffLkhd= [delta mean(prob1and2new,1)'./max(mean(prob1and2new,1))];  
