comparison Copy_of_multithreshold 1.46/bestFitPsychometicFunctions.m @ 28:02aa9826efe0

mainly multiThreshold

+% ------------------------------------------------------- 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
+