% 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