% Minimum Expected Entropy Staircase % % The staircase gives suggestions for which probe value to test next, % choosing the probe that will provide the most information (based on the % principle of minimum entropy = maximally unambiguous probability % distribution). Probes are chosen from a set of possible probe values % provided at staircase init, and their use is evaluated based on the % expected amount of information gain given a space of PSE and slope values % to test over. % % See MinExpEntStairDemo for an example and the comments in the method % functions below for use of the different staircase methods. % % By default, a psychometric function ranging from 0% to 100% is used, as % is suitable for discrimination experiments with a standard in the middle % of the possible stimulus parameter range. For other paradigms, such as % n-AFC detection tasks, one can set the guessrate input during staircase % init to 1/num_alternatives, e.g. .5 when doing a 2IFC detection task. % This guess rate is thus not the rate at which participants guess instead % of do your task (thats the lapse rate), it the minimum rate of correct % responses as determined by your design. NB: below discussion is based on % the default psychometric function with the full range, but all points are % equally valid for a scaled psychometric function. % % It is recommended to have the staircase determine the optimal next probe % based on only a random subset of the response history (see the % 'toggle_use_resp_subset' and 'toggle_use_resp_subset_prop' methods). This % makes its operation more robust for response errors and also avoids probe % oscillations when the fit estimate is converging. % When we are close to convergence, probes will tend to be near the 25% and % 75% points. If a probe is 25% and you answer '1' (pedestal faster, which % is likely, because it's near the correct 25% point), then for the next % trial the peak in expected entropy reduction will generally be the 75% % point, and vice versa. This can lead to undesirable probe sequences where % the correct response alternates 0,1,0,1,0,1. If you choose a random % subset, this will completely eliminate the problem. If the staircase has % converged to where there are two almost equal expected entropy minima, % then small variations due to the selection of subsets will randomly vary % which minimum emerges as lowest. % This strategy does not significantly affect optimal operation of the % staircase. Lots of probe values provide useful information. Therefore, it % is not crucial to have a highly accurate estimate of likelihoods, so % relatively few trials are sufficient (less than are needed to for final % estimates of PSE and DL). Throwing out trials for the staircase % computation yields robustness without much cost. % % Another option would be to load a non-uniform prior on the space of % possible location/mean/PSE and dispersion/slope parameters (known as mu % and sigma respectively for a cumulative Gaussian - see the loadprior % method). Probe sampling will then stay reasonable in early trials even if % there were a couple bad responses. But this strategy is not as robust as % using a random subset -- bad trials will continue to have an effect % throughout. % % In absence of anything to base the optimal probe value on, the first % probe is chosen randomly from the set of possible probes. When a prior % was loaded, a likelihood distribution is available based on which the % optional probe value can be computed. If for any other reason choosing % the next probe based on the measure of minimum expected entropy fails, % the staircase will fall back on the same random probe sampling strategy. % There is an option to set the first probe value to be tested (see the % set_first_value method), which, for the first trial only, will overrule % both of the above probe choice strategies. This can be useful if you want % to be sure that the first trial is an easy one so the participant knows % what to expect. % % Another measure for robustness is to choose a small lapse rate. If lapse % would be zero and a response error is made by the observer, immediately a % whole range of mean-slope combinations becomes impossible. If lapse rate % is non-zero, these would still have a non-zero probability and the % staircase can rebound. Therefore a lapse rate of 5% or even more % depending on task difficulty is always recommended. NB: in the default % discrimination setup of the staircase (guessrate is not specified or set % to 0), half of the lapse rate is taken off the bottom of the psychometric % function and half is taken off the top. So if the lapse rate is 0.05, the % psychometric function will range from 0.025 to 0.975. In the setup for a % n-AFC detection experiment when the psychometric function has a lower % bound of 1/num_alternatives, the lapse rate is taken off the top. So when % the guess_rate is set to .5 (2AFC) and the pase rate is set to .05, the % psychometric function will range from 0.05 to 0.95. % Note that the staircase does not support a 0 lapse rate in the first % place as it works with log-probability and we get in trouble if we would % take the log of a 0 probability. Any lapse rate lower than 1e-10 will be % adjusted to 1e-10 upon calling the init method. % % If the staircase gets stuck at one of the bounds of the probe set, check % that the sign of the slope space matches the expected sign of the % response. E.g., lets look at an experiment in which you are doing 2IFC % task in which the observer is asked to report which interval contained % the faster motion. If the observer choses the test over the pedestal % interval the response is 1, if the observer chosen the pedestal to be % faster, the response is 0. All slopes in the set would in this case be % positive as the low end of the probe space (slow speeds) is associated % with response 0 and the high end with response 1. If we however asked the % observer to indicate the slower interval, the slopes in our slope set % would not match the task, and the staircase would get stuck at one of the % probe bounds. In this case, the lower end of the probe space is % associated with the response 1 and the higher end with the response % 0--we'd thus have a negative slope for the fitted cumulative probability % function. % % The staircase currently only supports logistic and cumulative Gaussian % (default) psychometric functions (see the set_psychometric_func method), % but others could easily be implemented. Changes should be needed only to % the private fit_a_point method near the bottom of this mfile, providing % that the function is characterized by two parameters (which do not % necessarily have to be PSE and slope, though that is the terminology % here. % Should you implement such a function, please do send me your code to % dcnieho @at@ gmail.com. % % The above discussion assumes that response inputs to the process_resp % method are either 0 or 1 (see note above about their meaning) though in % practice anything larger than 0 is treated as 1 and anything lower than % 0, including 0, is treated as 0. the staircase can thus easily be % integrated with programs that use a 1, -1 response scheme. % % For actual offline fitting of your data, you would probably want to use a % dedicated toolbox such as Prins, N & Kingdom, F. A. A. (2009) Palamedes: % Matlab routines for analyzing psychophysical data. % http://www.palamedestoolbox.org. instead of using the function parameters % or PSE and DL returned from staircase methods get_fit and get_PSE_DL. % Also note that while the staircase runs far more robust when a small % lapse rate is assumed, it is common to either fit the psychometric % function without a lapse rate, or otherwise with the lapse rate as a free % parameter (possibily varying only over subjects, but not over conditions % within each subject). % % % References: % Based on the Minimum expected entropy staircase method developed by: % Saunders JA & Backus BT (2006). Perception of surface slant from % oriented textures. Journal of Vision 6(9), article 3 % % Discussions of conceptually similar staircases can be found in: % Kontsevich LL & Tyler CW (1999). Bayesian adaptive estimation of % psychometric slope and threshold. Vision Res 39(16), pp. 2729-37 % Lesmes LA, Lu ZL, Baek J & Albright TD (2010). Bayesian adaptive % estimation of the contrast sensitivity function: The quick CSF method. % Journal of Vision 10(3), article 17 % Copyright (c) 2011 by DC Niehorster and JA Saunders classdef MinExpEntStair < handle properties (Access=private, Hidden) % private member variables probeset = []; % possible probe values to be tested aset = []; % pse's tested (and fitted) bset = []; % slopes fitted agrid = []; bgrid = []; lapse_rate = []; % lapse/mistake rate guess_rate = []; % guess rate phist = []; % probe history rhist = double([]); % response history (0 or 1) loglik = []; lik = []; g0 = []; g1 = []; % likelihood lookup table qUseLookup = []; % can explicitly be set to true or false by user with likLookup = []; qLookupCompressed = false; % lots of overlap between likelihoods for different probe values, compute and store in a format making use of this overlap % option: use a subset of all data for choosing the next probe, default values: quse_subset = false; % use limited subset for computing next probe? Limited subset by discarding a fixed number of trials quse_subset_perc = false; % same as above, but instead use a percentage of the available data minsetsize = 10; % minimum size to start subsetting subsetsize = 3; % subset contains subsetsize less datapoints than full dataset percsetsize = .8; % percentage of data in set used % option: set the value to test if probe history is empty first_value = []; % first value to test instead of random or by prior % psychometric function that is used (default) psychofunc = []; psychofuncStr = 'cumGauss'; end % % subfunction % if nargin<1 || strcmpi(mode,'legacy') % fhndl = @MinExpEntStair_internal; % external_funs = {@init, @loadhistory, @loadprior, @toggle_use_resp_subset, @toggle_use_resp_subset_prop, @set_first_value, @set_use_lookup_table, @get_use_lookup_table, @set_psychometric_func, @get_psychometric_func, @get_next_probe, @process_resp, @get_history, @get_fit, @get_PSE_DL}; % external_funs_str = cellfun(@(x) strrep(func2str(x),[mfilename '/'],''),external_funs,'uni',false); % elseif strcmpi(mode,'v2') % % setup function handles % fhndl.init = @init; % fhndl.loadhistory = @loadhistory; % fhndl.loadprior = @loadprior; % fhndl.toggle_use_resp_subset = @toggle_use_resp_subset; % fhndl.toggle_use_resp_subset_prop = @toggle_use_resp_subset_prop; % fhndl.set_first_value = @set_first_value; % fhndl.set_use_lookup_table = @set_use_lookup_table; % fhndl.get_use_lookup_table = @get_use_lookup_table; % fhndl.set_psychometric_func = @set_psychometric_func; % fhndl.get_psychometric_func = @get_psychometric_func; % fhndl.get_next_probe = @get_next_probe; % fhndl.process_resp = @process_resp; % fhndl.get_history = @get_history; % fhndl.get_fit = @get_fit; % fhndl.get_PSE_DL = @get_PSE_DL; % % end % public interface methods function this = MinExpEntStair(mode) assert(nargin<1 || strcmpi(mode,'v2'),'Only mode v2 is still supported, ''legacy'' is no longer supported. You can just as well not provide this input.'); end % init function [] = init(this,probeset,meanset,slopeset,lapse_rate,guess_rate) this.probeset = probeset; this.aset = meanset; this.bset = slopeset; [this.agrid,this.bgrid] = meshgrid(this.aset,this.bset); % init with uniform probability, normalized this.loglik = zeros(size(this.agrid)) - log(numel(this.agrid)); this.lik = ones(size(this.agrid))./numel(this.agrid); % lapse rate and guess rate this.lapse_rate = lapse_rate; % the lapse rate cannot be exactly 0 as the computed % probability must not be exactly 0 so we can work with % log(prob) without trouble, so set it to 1e-10 at least. this.lapse_rate = max(this.lapse_rate,1e-10); % guess rate is optional, if not specified we assume a 2IFC % discrimination experiment where the guess rate is % irrelevant as function goes from always one option at the % one end to always the other option at the other end. if nargin<6 this.guess_rate = 0; else this.guess_rate = guess_rate; end % lapse rate: % 1. for a discrimination setup (guess_rate==0) the % lapserate basically means that instead of ranging from 0 % to 1, the psychometric function ranges from lapse_rate/2 % to 1-lapse_rate/2 % 2. for a detection setup, the lower bound is guess_rate % and the upper bound is 1-lapse_rate % lower bound of pyschometric function % and % range of pyschometric function if this.guess_rate==0 this.g0 = this.lapse_rate/2; this.g1 = 1 - this.lapse_rate; else this.g0 = this.guess_rate; this.g1 = 1 - this.lapse_rate - this.guess_rate; end this.set_psychometric_func(this.psychofuncStr); % calls precomputeLikelihoods() end %%% load bunch of previously run trials (need probes and %%% responses) function [] = loadhistory(this,probes,responses) this.phist = probes; this.rhist = responses; % refit likelihood up to this point this.fit_all(this.phist,this.rhist); end %%% load a prior likelihood, so that first probe is not chosen %%% randomly and you can influence evolution of the fit function [] = loadprior(this,priorlik) assert(all(this.loglik(:)==-log(numel(this.agrid))),'Cannot load prior if we have a likelihood already'); % this tests if it is not default inited assert(size(priorlik,1)==length(this.bset),'Number of rows in prior much match length of slope set') assert(size(priorlik,2)==length(this.aset),'Number of columns in prior much match length of mean set') assert(all(priorlik(:)>=0),'Loaded prior is not expected to be a log likelihood (that is: all your probabilities should be larger than or equal to 0!)'); this.loglik = normalize_loglik(log(priorlik)); this.lik = exp(this.loglik); end %%% use subset of data for computing next probe function [varargout] = toggle_use_resp_subset(this,minsetsize,subsetsize) % option: extract a probe and response subset for choosing % the next probe, and fit just those % when lots of trials ran, entropy function often has two % local minima, with their relative values switch back and % forth. This will lead to large oscillations in the probe % value being tested (one trial a probe from the beginning % of set, next trial a probe from the end and the from % beginning of set again). % We want to avoid these oscillations in probe values, % therefore we select a limited subset of data to calculate % the best next probe. this.quse_subset = ~this.quse_subset; assert(~(this.quse_subset && this.quse_subset_perc)); if nargin>1 % change defaults this.minsetsize = minsetsize; this.subsetsize = subsetsize; end varargout{1} = this.quse_subset; varargout{2} = this.minsetsize; varargout{3} = this.subsetsize; end %%% use subset of data for computing next probe function [varargout] = toggle_use_resp_subset_prop(this,minsetsize,percsetsize) % same as above, but now always use a proportion of the % available data this.quse_subset_perc = ~this.quse_subset_perc; assert(~(this.quse_subset_perc && this.quse_subset)); if nargin>1 % change defaults this.minsetsize = minsetsize; this.percsetsize = percsetsize; end varargout{1} = this.quse_subset_perc; varargout{2} = this.minsetsize; varargout{3} = this.percsetsize; end %%% set the first value to test. Normally the first is chosen %%% randomly or by using the prior that you loaded. If you prefer %%% to start at a fixed value, use this. function [] = set_first_value(this,first_value) this.first_value = first_value; if ~isempty(this.phist) warning('the first trial has already been run. Setting the first value now is pointless and it''ll be ignored'); end end %%% if set to true or false, for (not) using of a precomputed lookup %%% table instead of evaluating the psychometric function all the time. %%% call this before calling init as lookup table computation is %%% triggered at end of init function [] = set_use_lookup_table(this,qUseLookup) this.qUseLookup = qUseLookup; if this.qUseLookup && isempty(this.likLookup) this.precomputeLikelihoods(); end end %%% get if lookup table is currently used. function varargout = get_use_lookup_table(this) varargout{1} = this.qUseLookup; end %%% set the psychometric function to be used (default cumulative %%% Gaussian). Can be called at any time (but it will refit all %%% the data already present and thus remove the effect of any %%% priors). function [] = set_psychometric_func(this,funcID) if nargin>1 % currently supported: % 'cumGauss' - Cumulative Gaussian % 'logistic' - logistic function switch funcID case 'cumGauss' this.psychofunc = @(x,a,b) normcdf((x-a)./b); % 1 [ x - a ] % P = --- [ 1 + erf( ----------- ) ], % 2 [ b*sqrt(2) ] % where a and b are known as the mean (mu) and the standard % deviation (sigma) % http://en.wikipedia.org/wiki/Normal_distribution case 'logistic' this.psychofunc = @(x,a,b) 1./(1+exp(-(x-a)./b)); % 1 % P = ------------------, % -(x - a)/b % 1 + e^ % % where a and b are known as the mean (mu) and b is % proportional to the standard deviation (s) % http://en.wikipedia.org/wiki/Logistic_distribution otherwise error('Psychometric function "%s" not supported',funcID); end this.psychofuncStr = funcID; end % recompute lookup table this.precomputeLikelihoods(); % if there's any data already, refit it using the new % psychometric func. This would remove the effect of any % priors! if ~isempty(this.phist) ndata = min(length(this.phist),length(this.rhist)); this.fit_all(this.phist(1:ndata),this.rhist(1:ndata)); end end %%% get the psychometric function that is currently used. function [varargout] = get_psychometric_func(this) % currently possible outputs: % 'cumGauss' - Cumulative Gaussian % 'logistic' - logistic function varargout{1} = this.psychofuncStr; end %%% given history, get which probe is best to test next function [p,entexp,indmin] = get_next_probe(this) if isempty(this.phist) && ~isempty(this.first_value) % first trial and user requested a specific probe value to be tested p = this.first_value; [entexp,indmin] = deal([]); else [p,entexp,indmin] = this.getnextprobe(); if isempty(p) || isscalar(unique(this.loglik)) % if we couldn't compute expected entropy, or we have a % uniform likelihood on which calculation was based % (useless prior info, such as default inited), fall % back on random probe selection p = this.probeset(round(RandLim(1,1,length(this.probeset)))); [entexp,indmin] = deal([]); end end this.phist = [this.phist p]; end %%% fit likelihoods for new response function [] = process_resp(this,resp) % resp on current trial this.rhist(end+1) = resp; [this.loglik,this.lik] = this.fit_additional_data_point(this.loglik,this.phist(end),this.rhist(end)); end %%% retrieve probe and response history function [varargout] = get_history(this) varargout{1} = this.phist; varargout{2} = this.rhist; end %%% get fitted a (PSE) and b (slope) parameters and loglik. %%% This returns the fit of all data, also when subsetting is %%% enabled. function [varargout] = get_fit(this) kmin = find(this.loglik == max(this.loglik(:))); % most likely combination(s) of PSE and Slope varargout{1} = mean(this.agrid(kmin)); varargout{2} = mean(this.bgrid(kmin)); varargout{3} = this.loglik; end %%% get fitted PSE and DL (distance of 75% point from the 50% %%% point) and loglik. This returns the fit of all data, also %%% when subsetting is enabled. %%% This function is meant to be used for discrimination %%% experiments only (hence the terminology), although it will %%% return the inflection point and the distance between the %%% points that are equivalent to the 50% and 75% points after %%% scaling the psychometric function for all setups. function [varargout] = get_PSE_DL(this) [varargout{1:3}] = this.get_fit(); % convert b (dispersion) parameter to DL switch this.psychofuncStr case 'cumGauss' varargout{2} = varargout{2} * erfinv(.5)*sqrt(2); case 'logistic' varargout{2} = varargout{2} * log(3); otherwise error('Psychometric function "%s" not supported',this.psychofuncStr); end end end methods (Access=private, Hidden) % helpers (private functions, can only be called from the public % functions above) function [p,entexp,indmin] = getnextprobe(this) if length(this.rhist)>this.minsetsize && (this.quse_subset || this.quse_subset_perc) % select subset and fit if this.quse_subset_perc ind = NRandPerm(length(this.rhist),round(length(this.rhist)*this.percsetsize)); % select percentage of set else ind = NRandPerm(length(this.rhist),length(this.rhist)-this.subsetsize); % select set minus a few data points end [thellik,thelik] = this.fit_all(this.phist(ind),this.rhist(ind)); else % use likelihoods already fitted for all available data thelik = this.lik; thellik = this.loglik; end entexp = zeros(1,length(this.probeset)); for ksamp = 1:length(this.probeset) % p values for each possible model % these are used in multiple steps pvalsamp = this.fit_a_point(this.probeset(ksamp),1); % expected value is sum, weighted by lik pval = sum(pvalsamp(:).*thelik(:)); % two possibilities for next response, 0 or 1 % each would make a diff new likelihood function newloglik0 = thellik(:) + log(1 - pvalsamp(:)); newloglik1 = thellik(:) + log( pvalsamp(:)); % important! need to normalize newloglik0 = normalize_loglik(newloglik0); newloglik1 = normalize_loglik(newloglik1); % 0 and 1 for next response each has an entropy ent0 = sum(-exp(newloglik0).*newloglik0); ent1 = sum(-exp(newloglik1).*newloglik1); % probability pval of 0, probability (1-pval) of 1 % use these to get expected value of entropy entexp(ksamp) = ent0*(1-pval) + ent1*pval; end indmin = find(entexp == min(entexp),1); p = this.probeset(indmin); end function [loglik,lik] = fit_additional_data_point(this,loglik,probe,resp) % get likelihood of current point currlik = this.fit_a_point(probe,resp); % multiply with previous likelihoods loglik = loglik + log(currlik); % normalize loglik = normalize_loglik(loglik); if nargout>1 lik = exp(loglik); end end function [loglik,lik] = fit_all(this,probes,resps) if length(probes) ~= length(resps) error('Number of probe values and responses does not match'); end if strcmp(this.psychofuncStr,'cumGauss') % we have a fast one for this! loglik = FitCumGauss_MES(probes,resps,this.aset,this.bset,this.lapse_rate,this.guess_rate); % normalize loglik = normalize_loglik(loglik); else loglik = zeros(size(this.agrid)); for p=1:length(probes) loglik = this.fit_additional_data_point(loglik,probes(p),resps(p)); end % already normalized end lik = exp(loglik); end function pval = fit_a_point(this,probe,resp) if this.qUseLookup qProbe = this.probeset==probe; if this.qLookupCompressed pval = this.likLookup(:,[(end-length(this.aset)+1):end]-find(qProbe)+1); else pval = this.likLookup(:,:,qProbe); end else pval = this.evalLikelihood(probe); end % if response was "wrong", flip probs if resp <= 0 pval = 1-pval; end end function [] = precomputeLikelihoods(this) if isempty(this.aset) % called before init, parameter space not known yet, nothing to % do here return; end if ~isempty(this.qUseLookup) && ~this.qUseLookup % were not using lookup tables by users request, return return; end % determine if we want to precompute % first see if compressed format is possible. It is if same % stepsize for probeset and aset, as there is then significant % overlap between the pvalues for each probe level (could extend % this to one being multiples of the other...) stepP = mean(diff(this.probeset)); stepA = mean(diff(this.aset)); this.qLookupCompressed = abs(stepP-stepA)<=2*eps; % use lookup if compressed possible, or if table would be small, % or if user asked for it. if (isempty(this.qUseLookup) && (... this.qLookupCompressed || ... % same stepsize for probeset and aset numel(this.agrid)*length(this.probeset)/128/1024<3)... % small lookup table (by some arbitrary standard of what is small, which in this case is less than 3 mb) ) ||... (~isempty(this.qUseLookup) && this.qUseLookup) % user asked for it this.qUseLookup = true; nProbe = length(this.probeset); if this.qLookupCompressed [tempAGrid,tempBGrid] = meshgrid(linspace(this.probeset(1)-this.aset(1,end),this.probeset(end)-this.aset(1,1),length(this.aset)+length(this.probeset)-1),this.bset); this.likLookup = this.g0 + this.g1*this.psychofunc(0,tempAGrid,tempBGrid); else this.likLookup = zeros([size(this.agrid) nProbe]); for p=1:nProbe this.likLookup(:,:,p) = this.evalLikelihood(this.probeset(p)); end end else this.qUseLookup = false; end end function pval = evalLikelihood(this,probe) % evaluate psychometric function, incorporate lapse rate and guess rate pval = this.g0 + this.g1*this.psychofunc(probe,this.agrid,this.bgrid); end end end %% helpers function loglik = normalize_loglik(loglik) loglik = loglik - log(sum(exp(loglik(:)))); end