view multithreshold 1.46/bestFitPsychometicFunctions.m @ 38:c2204b18f4a2 tip

End nov big change

% ------------------------------------------------------- bestFitPsychometicFunctions
function [psy, levelsPhaseTwoBinVector, logistic, rareEvent]= bestFitPsychometicFunctions (levelsPhaseTwo, responsesPhaseTwo)
% bestFitPsychometicFunctions computes a psychometric function from a matrix of levels and responses
%output
%  psy is the empirical probabbility of a detection associated with levels in levelsPhaseTwoBinVector
%  logistic is a structure with results for fitting the logistic function
%    the logistic fit depends on whether maximum likelihood is used or least squares
%  rareEvent is a structure with results for fitting the rareEvent function
%   this is always calculated and lest squares is always used

global experiment stimulusParameters binFrequencies

% Generate a psychometic function by binning the levelsPhaseTwo
% x is a vector of [levelsPhaseTwo; responsesPhaseTwo]
% experiment.psyBinWidth defines the width of the bin
[psy, levelsPhaseTwoBinVector, binFrequencies,yesResponses, noResponses]= ...
    psychometricFunction(levelsPhaseTwo, responsesPhaseTwo, experiment.psyBinWidth);
% [levelsPhaseTwoBinVector; binFrequencies; psy];

% undefined slope - return
if isempty(psy)
    % slope is undefined
    psy=[]; levelsPhaseTwoBinVector=[];
    logistic.bestThreshold=NaN; logistic.bestK=NaN;  logistic.predictionsLOG=[]; logistic.predictionLevels=[];
    rareEvent.bestGain=NaN; rareEvent.bestVMin=NaN; rareEvent.predictionsRE=[]; rareEvent.predictionLevels=[];
    rareEvent.thresholddB= NaN; rareEvent.bestPaMindB=NaN;
    return
end

% find best fit rare event function
switch experiment.paradigm
    case 'gapDuration'
        % i.e. don't attempt this but visit the function to set null results
        rareEvent=fitRareEvent([], responsesPhaseTwo, stimulusParameters.targetDuration);
    otherwise
        rareEvent=fitRareEvent(levelsPhaseTwo, responsesPhaseTwo, stimulusParameters.targetDuration);
end

% find best logistic fit
logisticFunc=' 1./(1+exp(-a2.*(x-a1)));';
switch experiment.functionEstMethod
    % least squares estimate
    case {'logisticLS', 'rareEvent','peaksAndTroughs'}
        [a1, a2, Euclid]=fitFunctionUsingLeastSquares(levelsPhaseTwo, responsesPhaseTwo, ...
            logisticFunc, experiment.possLogSlopes, experiment.meanSearchStep);
%         [a1, a2, Euclid]=fitLogisticUsingLS(psy, levelsPhaseTwoBinVector); %! does not give same result
        logistic.bestThreshold=a1;
        logistic.bestK=a2;
        logistic.predictionsLOG= ...
            1./(1+exp(-logistic.bestK.*(experiment.predictionLevels-logistic.bestThreshold)));
        logistic.predictionLevels=experiment.predictionLevels;
        logistic.Euclid = Euclid;
        %         predResponses=1./(1+exp(-logistic.bestK.*(levelsPhaseTwo-logistic.bestThreshold)));
        %         [levelsPhaseTwo' responsesPhaseTwo' predResponses' responsesPhaseTwo'-predResponses' ]

        % maximum likelihood fitting
    case 'logisticML'
        [a1, a2, Euclid]=fitFunctionUsingMaxLikelihood (levelsPhaseTwoBinVector,...
            psy, yesResponses, noResponses, logisticFunc, experiment.possLogSlopes,experiment.meanSearchStep);
        logistic.bestThreshold=a1;
        logistic.bestK=a2;
        logistic.predictionsLOG= ...
            1./(1+exp(-logistic.bestK.*(experiment.predictionLevels-logistic.bestThreshold)));
        logistic.predictionLevels=experiment.predictionLevels;
        logistic.Euclid = Euclid;
end

% disp(num2str([logistic.bestThreshold logistic.bestK logistic.Euclid]))
% disp(num2str([ rareEvent.thresholddB rareEvent.bestGain rareEvent.bestVMin rareEvent.Euclid]))

% --------------------------------------------- fitLogisticUsingLS
function [mean, slope, Euclid]=fitLogisticUsingLS(p,L)

%!! each bin needs to be weighted by its size 

idx1=find(p>0);
p=p(idx1);
L=L(idx1);
idx1=find(p<1);
p=p(idx1);
L=L(idx1);

if length(L)<2
    mean=NaN;
    slope=NaN;
    Euclid=NaN;
    return
end

y=L;
x=log(1./p -1);
a =polyfit(x, y, 1);
mean=a(2);
slope=-1/a(1);

% euclid
predy=polyval(a,x);
Euclid=sum((predy-y).^2);
[mean slope Euclid];

% --------------------------------------------- fitFunctionUsingLeastSquares
function [a1, a2, Euclid]=fitFunctionUsingLeastSquares(levelsPhaseTwo, responsesPhaseTwo, func, possSlopes, meanSearchStep)

% nResponses= yesResponses+noResponses;
% idx=find(not(nResponses==0));
% levelsPhaseTwo=levelsPhaseTwo(idx);
% yesResponses=yesResponses(idx);
% noResponses=noResponses(idx);
% nResponses=nResponses(idx);

% psy=yesResponses./nResponses;

possThresholds=min(levelsPhaseTwo):meanSearchStep:max(levelsPhaseTwo);


% create vectors representing all combinations of a1 and a2
nPossThresh=length(possThresholds);
nPossSlopes=length(possSlopes);
a1=repmat(possThresholds',nPossSlopes,1);
a2=repmat(possSlopes,nPossThresh,1);
a2=reshape(a2,nPossThresh*nPossSlopes,1);

% each response is predicted by every possible threshold and slope
LS=ones(size(a1));
for i=1:length(levelsPhaseTwo)
    x=levelsPhaseTwo(i);
    eval(['y=' func]);
    actual=responsesPhaseTwo(i);
    error= (actual - y).^2;
    LS=LS+error;
end

[Euclid, idx]=min(LS);

a1=a1(idx);
a2=a2(idx);

% x=levelsPhaseTwo;
% eval(['y=' func]);
% Euclid= sum((psy-y).^2);

% plot(levelsPhaseTwo,y,'r')
% hold on
% plot(levelsPhaseTwo,psy,'o')

% --------------------------------------------- fitFunctionUsingMaxLikelihood
function [a1, a2, Euclid]=fitFunctionUsingMaxLikelihood...
    (levelsPhaseTwo,psy, yesResponses, noResponses, func, possible_a2, meanSearchStep)

% fitFunctionUsingMaxLikelihood fits the function in 'func' to binned yes-no data.
% levelsPhaseTwo specifies the bin centers
% yesResponses and noResponses are the niumber of yes and no responsesPhaseTwo in each bin
% 'func' is a function of a1, b1, x; e.g. func='1./(1+exp(-a2.*(x-a1)));';
%  and a2, a1 are parameters to be discovered.
% If bins are empty (i.e. neither yes or no), these are eliminated
%

nResponses=yesResponses+noResponses;
possible_a1=min(levelsPhaseTwo):meanSearchStep:max(levelsPhaseTwo);

% if nargin<6
%     possible_a2=-2:.05:4;
% end

% create vectors representing all combinations of a1 and a2
nPossThresh=length(possible_a1);
nPossSlopes=length(possible_a2);
a1=repmat(possible_a1',nPossSlopes,1);
a2=repmat(possible_a2,nPossThresh,1);
a2=reshape(a2,nPossThresh*nPossSlopes,1);

cumulativeProbability=ones(1, length(a2));
cumulativeProbability=cumulativeProbability';

for i=1:length(levelsPhaseTwo)
    x=levelsPhaseTwo(i);
    eval(['y=' func]);
    funcVal=(1-y).^yesResponses(i);
    cumulativeProbability= cumulativeProbability.*funcVal;
    funcVal=y.^noResponses(i);
    cumulativeProbability= cumulativeProbability.*funcVal;
end

[maxProb idx]=max(cumulativeProbability);
a1=a1(idx);
a2=a2(idx);

% x=levelsPhaseTwo;
eval(['y=' func]);
Euclid= sum(nResponses.*(psy-y).^2);

% figure(1), clf
% plot(levelsPhaseTwo,y)
% hold on
% plot(levelsPhaseTwo,psy,'o')



% --------------------------------------------------- fitRareEvent
function rareEvent=fitRareEvent(stimulusLevels, responsesPhaseTwo, duration, gains, Vmins)
% least squares estimate of *rare event* function
% model is: r = g P – A   ... g=events/s/Pa, A= events/s, P= pressure (Pa)
% psy=1- exp(d(gP –A  ))   ... d=duration
% p(event)=gain*levelmPa -Vmin
% 'responsesPhaseTwo' is a binary vector of subject's decision.
% 'stimulusLevels' are the corresponding signal levesl (values)
% duration is required to compute the expectation of an event occurring
% gains is an optional list of gains to be tried
% Vmins is an optional list of Vmins to be tried

global experiment
if nargin<5
    minVmin=.1; maxVmin=10; nVmins=100;
    Vmins=[0 logspace(log10(minVmin),log10(maxVmin),nVmins)];
    % disp(Vmins)
    gainMin=0.0001; gainMax=1; nGains=100; 
    gains=[1e-5 logspace(log10(gainMin),log10(gainMax),nGains)];
end

rareEvent.bestGain=NaN;
rareEvent.bestVMin=NaN;
rareEvent.thresholddB=0;
rareEvent.bestPaMindB=NaN;
rareEvent.predictionLevels=[];
rareEvent.predictionsRE=[];
rareEvent.Euclid=NaN;

if isempty(stimulusLevels), return, end

% expected slope is negative, rareEvent can not be used
if experiment.psyFunSlope<0
    return
end

% NB calculations in microPascals!
stimulusLevelsAsPressure=28 * 10.^(stimulusLevels/20);

% allGains=reshape(repmat(gains,nVmins,1), 1, nVmins*nGains);
% allVmins=repmat(Vmins, 1, nGains);

    predictions=NaN*zeros(1,length(stimulusLevels));
    gainCount=0;
    Euclid=inf; bestVmin=0; bestGain=0;
    for gain= gains
        gainCount=gainCount+1;
        VminCount=0;
        
        % 1 – exp(-d (g P – A))
        for Vmin=Vmins
            VminCount=VminCount+1;
            % all levelsPhaseTwo are simultaneously assessed
            
            % original linear function
%             gP_Vmin=gain*stimulusLevelsAsPressure-Vmin;
%             idx=(gP_Vmin>0);
            %   predictions(idx)= 1-exp(-duration*(gP_Vmin(idx)));
            %   new log function log(r) =gP-A
            
            % log function
            gP_Vmin=gain*stimulusLevelsAsPressure-Vmin;
            idx=(gP_Vmin>0);
            predictions(idx)= 1-exp(-duration*exp(gP_Vmin(idx)));
            predictions(~idx)=0;
            
%             % square function
%             P_Vmin=stimulusLevelsAsPressure-Vmin;
%             idx=(P_Vmin)>0;
%             predictions(idx)= 1-exp(-duration*gain*(P_Vmin(idx)).^2);
%             predictions(~idx)=0;
            
            
            % NB the error is equally weighted for each stimulus/response pair
            % this is not the same as equal weighting for each distinct level
            error=(predictions - responsesPhaseTwo).^2;
            error=mean(error(~isnan(error)));
            if error<Euclid
                Euclid=error;
                bestVmin=Vmin;
                bestVminCount=VminCount;
                %                 bestGainCount=gainCount;
                bestGain=gain;
            end
        end
    end
    % disp(Vmins)

% if bestGainCount==1 | bestGainCount==nGains
%     disp(['gain estimate ' num2str(gains(bestGainCount)) ' may be out of range'])
% end

[rareEvent.Euclid idx]=min(Euclid);
rareEvent.bestGain=bestGain;
rareEvent.bestVMin=bestVmin;
rareEvent.thresholdPa=(-log(0.5)/duration + rareEvent.bestVMin)/rareEvent.bestGain;
rareEvent.thresholddB=20*log10(rareEvent.thresholdPa/28);
rareEvent.bestPaMindB=20*log10((rareEvent.bestVMin/rareEvent.bestGain)/28);

predictionLevels= -50:1:120;
rareEvent.predictionLevels=predictionLevels;
stimulusLevelsAsPressure=28*10.^(predictionLevels/20);

% gP_Vmin=rareEvent.bestGain*stimulusLevelsAsPressure-rareEvent.bestVMin;
% rareEvent.predictionsRE= 1-exp(-duration*(gP_Vmin));

gP_Vmin=rareEvent.bestGain*stimulusLevelsAsPressure-rareEvent.bestVMin;
rareEvent.predictionsRE=  1-exp(-duration*exp(gP_Vmin));

%             P_Vmin=stimulusLevelsAsPressure-Vmin;
%             predictions= 1-exp(-duration*gain*P_Vmin.^2);
% predictions(predictions<0)=0;
%             rareEvent.predictionsRE=predictions;

% rareEvent