rmeddis@28: % ------------------------------------------------------- bestFitPsychometicFunctions rmeddis@28: function [psy, levelsPhaseTwoBinVector, logistic, rareEvent]= bestFitPsychometicFunctions (levelsPhaseTwo, responsesPhaseTwo) rmeddis@28: % bestFitPsychometicFunctions computes a psychometric function from a matrix of levels and responses rmeddis@28: %output rmeddis@28: % psy is the empirical probabbility of a detection associated with levels in levelsPhaseTwoBinVector rmeddis@28: % logistic is a structure with results for fitting the logistic function rmeddis@28: % the logistic fit depends on whether maximum likelihood is used or least squares rmeddis@28: % rareEvent is a structure with results for fitting the rareEvent function rmeddis@28: % this is always calculated and lest squares is always used rmeddis@28: rmeddis@28: global experiment stimulusParameters binFrequencies rmeddis@28: rmeddis@28: % Generate a psychometic function by binning the levelsPhaseTwo rmeddis@28: % x is a vector of [levelsPhaseTwo; responsesPhaseTwo] rmeddis@28: % experiment.psyBinWidth defines the width of the bin rmeddis@28: [psy, levelsPhaseTwoBinVector, binFrequencies,yesResponses, noResponses]= ... rmeddis@28: psychometricFunction(levelsPhaseTwo, responsesPhaseTwo, experiment.psyBinWidth); rmeddis@28: % [levelsPhaseTwoBinVector; binFrequencies; psy]; rmeddis@28: rmeddis@28: % undefined slope - return rmeddis@28: if isempty(psy) rmeddis@28: % slope is undefined rmeddis@28: psy=[]; levelsPhaseTwoBinVector=[]; rmeddis@28: logistic.bestThreshold=NaN; logistic.bestK=NaN; logistic.predictionsLOG=[]; logistic.predictionLevels=[]; rmeddis@28: rareEvent.bestGain=NaN; rareEvent.bestVMin=NaN; rareEvent.predictionsRE=[]; rareEvent.predictionLevels=[]; rmeddis@28: rareEvent.thresholddB= NaN; rareEvent.bestPaMindB=NaN; rmeddis@28: return rmeddis@28: end rmeddis@28: rmeddis@28: % find best fit rare event function rmeddis@28: switch experiment.paradigm rmeddis@28: case 'gapDuration' rmeddis@28: % i.e. don't attempt this but visit the function to set null results rmeddis@28: rareEvent=fitRareEvent([], responsesPhaseTwo, stimulusParameters.targetDuration); rmeddis@28: otherwise rmeddis@28: rareEvent=fitRareEvent(levelsPhaseTwo, responsesPhaseTwo, stimulusParameters.targetDuration); rmeddis@28: end rmeddis@28: rmeddis@28: % find best logistic fit rmeddis@28: logisticFunc=' 1./(1+exp(-a2.*(x-a1)));'; rmeddis@28: switch experiment.functionEstMethod rmeddis@28: % least squares estimate rmeddis@28: case {'logisticLS', 'rareEvent','peaksAndTroughs'} rmeddis@28: [a1, a2, Euclid]=fitFunctionUsingLeastSquares(levelsPhaseTwo, responsesPhaseTwo, ... rmeddis@28: logisticFunc, experiment.possLogSlopes, experiment.meanSearchStep); rmeddis@28: % [a1, a2, Euclid]=fitLogisticUsingLS(psy, levelsPhaseTwoBinVector); %! does not give same result rmeddis@28: logistic.bestThreshold=a1; rmeddis@28: logistic.bestK=a2; rmeddis@28: logistic.predictionsLOG= ... rmeddis@28: 1./(1+exp(-logistic.bestK.*(experiment.predictionLevels-logistic.bestThreshold))); rmeddis@28: logistic.predictionLevels=experiment.predictionLevels; rmeddis@28: logistic.Euclid = Euclid; rmeddis@28: % predResponses=1./(1+exp(-logistic.bestK.*(levelsPhaseTwo-logistic.bestThreshold))); rmeddis@28: % [levelsPhaseTwo' responsesPhaseTwo' predResponses' responsesPhaseTwo'-predResponses' ] rmeddis@28: rmeddis@28: % maximum likelihood fitting rmeddis@28: case 'logisticML' rmeddis@28: [a1, a2, Euclid]=fitFunctionUsingMaxLikelihood (levelsPhaseTwoBinVector,... rmeddis@28: psy, yesResponses, noResponses, logisticFunc, experiment.possLogSlopes,experiment.meanSearchStep); rmeddis@28: logistic.bestThreshold=a1; rmeddis@28: logistic.bestK=a2; rmeddis@28: logistic.predictionsLOG= ... rmeddis@28: 1./(1+exp(-logistic.bestK.*(experiment.predictionLevels-logistic.bestThreshold))); rmeddis@28: logistic.predictionLevels=experiment.predictionLevels; rmeddis@28: logistic.Euclid = Euclid; rmeddis@28: end rmeddis@28: rmeddis@28: % disp(num2str([logistic.bestThreshold logistic.bestK logistic.Euclid])) rmeddis@28: % disp(num2str([ rareEvent.thresholddB rareEvent.bestGain rareEvent.bestVMin rareEvent.Euclid])) rmeddis@28: rmeddis@28: % --------------------------------------------- fitLogisticUsingLS rmeddis@28: function [mean, slope, Euclid]=fitLogisticUsingLS(p,L) rmeddis@28: rmeddis@28: %!! each bin needs to be weighted by its size rmeddis@28: rmeddis@28: idx1=find(p>0); rmeddis@28: p=p(idx1); rmeddis@28: L=L(idx1); rmeddis@28: idx1=find(p<1); rmeddis@28: p=p(idx1); rmeddis@28: L=L(idx1); rmeddis@28: rmeddis@28: if length(L)<2 rmeddis@28: mean=NaN; rmeddis@28: slope=NaN; rmeddis@28: Euclid=NaN; rmeddis@28: return rmeddis@28: end rmeddis@28: rmeddis@28: y=L; rmeddis@28: x=log(1./p -1); rmeddis@28: a =polyfit(x, y, 1); rmeddis@28: mean=a(2); rmeddis@28: slope=-1/a(1); rmeddis@28: rmeddis@28: % euclid rmeddis@28: predy=polyval(a,x); rmeddis@28: Euclid=sum((predy-y).^2); rmeddis@28: [mean slope Euclid]; rmeddis@28: rmeddis@28: % --------------------------------------------- fitFunctionUsingLeastSquares rmeddis@28: function [a1, a2, Euclid]=fitFunctionUsingLeastSquares(levelsPhaseTwo, responsesPhaseTwo, func, possSlopes, meanSearchStep) rmeddis@28: rmeddis@28: % nResponses= yesResponses+noResponses; rmeddis@28: % idx=find(not(nResponses==0)); rmeddis@28: % levelsPhaseTwo=levelsPhaseTwo(idx); rmeddis@28: % yesResponses=yesResponses(idx); rmeddis@28: % noResponses=noResponses(idx); rmeddis@28: % nResponses=nResponses(idx); rmeddis@28: rmeddis@28: % psy=yesResponses./nResponses; rmeddis@28: rmeddis@28: possThresholds=min(levelsPhaseTwo):meanSearchStep:max(levelsPhaseTwo); rmeddis@28: rmeddis@28: rmeddis@28: % create vectors representing all combinations of a1 and a2 rmeddis@28: nPossThresh=length(possThresholds); rmeddis@28: nPossSlopes=length(possSlopes); rmeddis@28: a1=repmat(possThresholds',nPossSlopes,1); rmeddis@28: a2=repmat(possSlopes,nPossThresh,1); rmeddis@28: a2=reshape(a2,nPossThresh*nPossSlopes,1); rmeddis@28: rmeddis@28: % each response is predicted by every possible threshold and slope rmeddis@28: LS=ones(size(a1)); rmeddis@28: for i=1:length(levelsPhaseTwo) rmeddis@28: x=levelsPhaseTwo(i); rmeddis@28: eval(['y=' func]); rmeddis@28: actual=responsesPhaseTwo(i); rmeddis@28: error= (actual - y).^2; rmeddis@28: LS=LS+error; rmeddis@28: end rmeddis@28: rmeddis@28: [Euclid, idx]=min(LS); rmeddis@28: rmeddis@28: a1=a1(idx); rmeddis@28: a2=a2(idx); rmeddis@28: rmeddis@28: % x=levelsPhaseTwo; rmeddis@28: % eval(['y=' func]); rmeddis@28: % Euclid= sum((psy-y).^2); rmeddis@28: rmeddis@28: % plot(levelsPhaseTwo,y,'r') rmeddis@28: % hold on rmeddis@28: % plot(levelsPhaseTwo,psy,'o') rmeddis@28: rmeddis@28: % --------------------------------------------- fitFunctionUsingMaxLikelihood rmeddis@28: function [a1, a2, Euclid]=fitFunctionUsingMaxLikelihood... rmeddis@28: (levelsPhaseTwo,psy, yesResponses, noResponses, func, possible_a2, meanSearchStep) rmeddis@28: rmeddis@28: % fitFunctionUsingMaxLikelihood fits the function in 'func' to binned yes-no data. rmeddis@28: % levelsPhaseTwo specifies the bin centers rmeddis@28: % yesResponses and noResponses are the niumber of yes and no responsesPhaseTwo in each bin rmeddis@28: % 'func' is a function of a1, b1, x; e.g. func='1./(1+exp(-a2.*(x-a1)));'; rmeddis@28: % and a2, a1 are parameters to be discovered. rmeddis@28: % If bins are empty (i.e. neither yes or no), these are eliminated rmeddis@28: % rmeddis@28: rmeddis@28: nResponses=yesResponses+noResponses; rmeddis@28: possible_a1=min(levelsPhaseTwo):meanSearchStep:max(levelsPhaseTwo); rmeddis@28: rmeddis@28: % if nargin<6 rmeddis@28: % possible_a2=-2:.05:4; rmeddis@28: % end rmeddis@28: rmeddis@28: % create vectors representing all combinations of a1 and a2 rmeddis@28: nPossThresh=length(possible_a1); rmeddis@28: nPossSlopes=length(possible_a2); rmeddis@28: a1=repmat(possible_a1',nPossSlopes,1); rmeddis@28: a2=repmat(possible_a2,nPossThresh,1); rmeddis@28: a2=reshape(a2,nPossThresh*nPossSlopes,1); rmeddis@28: rmeddis@28: cumulativeProbability=ones(1, length(a2)); rmeddis@28: cumulativeProbability=cumulativeProbability'; rmeddis@28: rmeddis@28: for i=1:length(levelsPhaseTwo) rmeddis@28: x=levelsPhaseTwo(i); rmeddis@28: eval(['y=' func]); rmeddis@28: funcVal=(1-y).^yesResponses(i); rmeddis@28: cumulativeProbability= cumulativeProbability.*funcVal; rmeddis@28: funcVal=y.^noResponses(i); rmeddis@28: cumulativeProbability= cumulativeProbability.*funcVal; rmeddis@28: end rmeddis@28: rmeddis@28: [maxProb idx]=max(cumulativeProbability); rmeddis@28: a1=a1(idx); rmeddis@28: a2=a2(idx); rmeddis@28: rmeddis@28: % x=levelsPhaseTwo; rmeddis@28: eval(['y=' func]); rmeddis@28: Euclid= sum(nResponses.*(psy-y).^2); rmeddis@28: rmeddis@28: % figure(1), clf rmeddis@28: % plot(levelsPhaseTwo,y) rmeddis@28: % hold on rmeddis@28: % plot(levelsPhaseTwo,psy,'o') rmeddis@28: rmeddis@28: rmeddis@28: rmeddis@28: % --------------------------------------------------- fitRareEvent rmeddis@28: function rareEvent=fitRareEvent(stimulusLevels, responsesPhaseTwo, duration, gains, Vmins) rmeddis@28: % least squares estimate of *rare event* function rmeddis@28: % model is: r = g P – A ... g=events/s/Pa, A= events/s, P= pressure (Pa) rmeddis@28: % psy=1- exp(d(gP –A )) ... d=duration rmeddis@28: % p(event)=gain*levelmPa -Vmin rmeddis@28: % 'responsesPhaseTwo' is a binary vector of subject's decision. rmeddis@28: % 'stimulusLevels' are the corresponding signal levesl (values) rmeddis@28: % duration is required to compute the expectation of an event occurring rmeddis@28: % gains is an optional list of gains to be tried rmeddis@28: % Vmins is an optional list of Vmins to be tried rmeddis@28: rmeddis@28: global experiment rmeddis@28: if nargin<5 rmeddis@28: minVmin=.1; maxVmin=10; nVmins=100; rmeddis@28: Vmins=[0 logspace(log10(minVmin),log10(maxVmin),nVmins)]; rmeddis@28: % disp(Vmins) rmeddis@28: gainMin=0.0001; gainMax=1; nGains=100; rmeddis@28: gains=[1e-5 logspace(log10(gainMin),log10(gainMax),nGains)]; rmeddis@28: end rmeddis@28: rmeddis@28: rareEvent.bestGain=NaN; rmeddis@28: rareEvent.bestVMin=NaN; rmeddis@28: rareEvent.thresholddB=0; rmeddis@28: rareEvent.bestPaMindB=NaN; rmeddis@28: rareEvent.predictionLevels=[]; rmeddis@28: rareEvent.predictionsRE=[]; rmeddis@28: rareEvent.Euclid=NaN; rmeddis@28: rmeddis@28: if isempty(stimulusLevels), return, end rmeddis@28: rmeddis@28: % expected slope is negative, rareEvent can not be used rmeddis@28: if experiment.psyFunSlope<0 rmeddis@28: return rmeddis@28: end rmeddis@28: rmeddis@28: % NB calculations in microPascals! rmeddis@28: stimulusLevelsAsPressure=28 * 10.^(stimulusLevels/20); rmeddis@28: rmeddis@28: % allGains=reshape(repmat(gains,nVmins,1), 1, nVmins*nGains); rmeddis@28: % allVmins=repmat(Vmins, 1, nGains); rmeddis@28: rmeddis@28: predictions=NaN*zeros(1,length(stimulusLevels)); rmeddis@28: gainCount=0; rmeddis@28: Euclid=inf; bestVmin=0; bestGain=0; rmeddis@28: for gain= gains rmeddis@28: gainCount=gainCount+1; rmeddis@28: VminCount=0; rmeddis@28: rmeddis@28: % 1 – exp(-d (g P – A)) rmeddis@28: for Vmin=Vmins rmeddis@28: VminCount=VminCount+1; rmeddis@28: % all levelsPhaseTwo are simultaneously assessed rmeddis@28: rmeddis@28: % original linear function rmeddis@28: % gP_Vmin=gain*stimulusLevelsAsPressure-Vmin; rmeddis@28: % idx=(gP_Vmin>0); rmeddis@28: % predictions(idx)= 1-exp(-duration*(gP_Vmin(idx))); rmeddis@28: % new log function log(r) =gP-A rmeddis@28: rmeddis@28: % log function rmeddis@28: gP_Vmin=gain*stimulusLevelsAsPressure-Vmin; rmeddis@28: idx=(gP_Vmin>0); rmeddis@28: predictions(idx)= 1-exp(-duration*exp(gP_Vmin(idx))); rmeddis@28: predictions(~idx)=0; rmeddis@28: rmeddis@28: % % square function rmeddis@28: % P_Vmin=stimulusLevelsAsPressure-Vmin; rmeddis@28: % idx=(P_Vmin)>0; rmeddis@28: % predictions(idx)= 1-exp(-duration*gain*(P_Vmin(idx)).^2); rmeddis@28: % predictions(~idx)=0; rmeddis@28: rmeddis@28: rmeddis@28: % NB the error is equally weighted for each stimulus/response pair rmeddis@28: % this is not the same as equal weighting for each distinct level rmeddis@28: error=(predictions - responsesPhaseTwo).^2; rmeddis@28: error=mean(error(~isnan(error))); rmeddis@28: if error