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