Fitting a better psychometric curve (Wichmann and Hill 2001)

Code
GitHub: Matlab code and Python package

Introduction
Psychometric curves quantify the relationship between stimulus parameters and an observer's subjective experience. For an introduction to the concept of psychometric curves in signal analysis theory and a function and instructions to fit a basic curve in Matlab, see Introduction to psychometric curves and fitting a simple psychometric function. This article will focus on fitting a more complex curve with additional parameters that allow for the subjects fallibility, which improves estimation of the key parameter, discrimination sensitivity.

The psychometric function used here is from the paper The psychometric function:I. Fitting, sampling, and goodness of fit, F. Wichmann and N. Hill 2001. This paper, and its companion paper, provide detailed theory and methods for fitting a psychometric function, assessing its goodness of fit, and finding confidence intervals for the ranges of the parameters estimated. This article will demonstrate how to fit Wichmann and Hill's psychometric function in Matlab using the same data from a fictional 2AFC experiment in this post. This will allow an interested reader to make a direct comparison between the Wichmann and Hill curve and the approach using Matlab's glmfit function with a built in binomial distribution and a logit link function.




The Wichmann and Hill psychometric function
The psychometric function presented in Wichmann and Hill's paper is a cumulative gaussian function with 4 parameters:

• Mean (u): The mean value of the distribution representing subject bias.
• Standard deviation (v): The variation of the distribution representing the subjects discrimination sensitivity.
• Guess rate (g) and lapse rate (l): Two additional parameters representing the subjects fallibility (ie. potential inability to ever reach 100% performance) at each end of the distribution/stimulus spectrum. 

The parameters g and l constrain the limits of the cumulative distribution that provides the sigmoid shape for the psychometric curve:

y = g+(1-g-l)*dist,

where dist is the distribution to fit, which is a cumulative gaussian

dist = 0.5*(1+erf((x-u)/sqrt(2*v^2)))

In Matlab we will define this function with the code

F=@(g,l,u,v,x) g+(1-g-l)*0.5*(1+erf((x-u)/sqrt(2*v^2)));

This creates an anonymous function in the workspace that expects the 4 parameters as inputs (g,l,u,v) along with an x axis (x). This function is then passed to Matlab's fit function, which will fit the curve to our data. For example,

ffit=fit(xAxis,yData,F)

The defining of the psychometric function (F) and a few other housekeeping tasks are handled by FitPsycheCurveWH.m.

Fitting the curve
We'll reuse the data and 2 alternative forced choice (2AFC) experiment from Introduction to psychometric curves and fitting a simple psychometric function, fit some curves using FitPsycheCurveWH.m, and calculate the bias and threshold under different conditions. Imagine we have the following conditions in which we want to test the subjects ability to discriminate between the colours white and black. The subject is presented with a stimulus and must respond by categorising it as either "black" or "white".

• Neural background - the "default" condition with no distractions.
• Dark background - the target stimulus is presented with a dark background that contrasts with white stimuli, perhaps making them look whiter. This increases the subjects likelihood of categorising a stimulus as white, and creating a bias (at least it does in our made up data anyway!).
• Obscured target - the stimulus is obscured in some way, making the task harder and increasing the subjects threshold (decreasing discrimination sensitivity); meaning a greater difference in shade would be required to make the white/black categorisation as accurately as in the default condition. 

% Make up some data

x=1:10
yNeutralBg = [5, 5, 10, 20, 40, 60, 80, 90, 95, 95]/100
yDarkBg = [0, 0, 2.5, 5, 10, 20, 40, 60, 80, 90]/100
yObscured = [25, 25, 30, 35, 45, 55, 65, 70, 75, 75]/100

% Plot the made up data
figure
scatter(x,yNeutralBg)
hold on
scatter(x,yDarkBg)
scatter(x,yObscured)
ylabel('Proportion "black" responses')
xlabel('Greyscale, 1=white, 10=black')

% Fit for neutral background
SPs = [0.35, 0.1, 10, 10; % Upper limits for g, l, u ,v
    0.01, 0.05, 5, 1; % Start points for g, l, u ,v
    0, 0, 0, 0]; % Lower limits for  g, l, u ,v

FitPsycheCurveWH takes a variable number of inputs; either just the x axis and y data, or the x axis, y data, and estimations for the 4 coefficients. The fitting itself is done by Matlab's fit function, which accepts start points and upper and lower limits for each of the coefficients. These values narrow down the parameter space that needs to be explored, and setting sensible values before fitting can reduce processing and increases the chances of a sensible fit being found. Here these values are entered in the 3x4 matrix SPs, where there's a column for each coefficient (g, l, u, v, respectively). Row 1 contains the upper limits, row 2 the start points, and row 3 the lower limits.

[coeffsNeutralBg, curveNeutralBg] = ...
    FitPsycheCurveWH(x, yNeutralBg, SPs);

Similar to the FitPsychCurveLogit function, FitPsycheCurveWH outputs the fitted coefficients (here in the variable coeffsNeutralBG) and a higher resolution curve for plotting (here in curveNeutralBg). The coefficients are contained in a fit object created by the fit function, which is evaluated by feval to get the data for the higher resolution curve. As a convenient alternative, fit objects can also be passed directly to functions like plot, where they will be evaluated over a specified or automatic x axis first, then plotted.

% Fit for dark background
SPs = [0.35, 0.35, inf, inf;
    0.01, 0.01, 7.5, 1;
    0, 0, 0, 0];
[coeffsyDarkBg, curveDarkBg] = ...
    FitPsycheCurveWH(x, yDarkBg, SPs);
% Fit for obscured target
SPs = [0.5, 0.5, inf, inf;
    0.01, 0.01, 7.5, 1;
    0, 0, 0, 0];
[coeffsObscured, curveObscured] = ...
    FitPsycheCurveWH(x, yObscured, SPs);

% Plot psychometic curves
plot(curveNeutralBg(:,1), curveNeutralBg(:,2), 'LineStyle', '--')
plot(curveDarkBg(:,1), curveDarkBg(:,2), 'LineStyle', '--')
plot(curveObscured(:,1), curveObscured(:,2), 'LineStyle', '--')
legend('Neutral background', 'Dark background', 'Obscured target', ...
'Neutral background fit', 'Dark background fit', 'Obscured target fit');

These look like nice fits. Admittedly this data is cleaner than most real data, but comparing the fits to those generate by FitPsychCurveLogit, there's noticably less deviation from the curves even with such clean data, and particularly in the case of the obscured target condition.

Bias/PSE and discrimination sensitivity
Conveniently, for this psychometric function, the coefficient values u (distribution mean) and v (distribution standard deviation) directly represent the bias and discrimination threshold, respectively, so no additional calculations are necessary. These values are available in the fit objects, and can be accessed using the .notion, for example:

coeffsyDarkBg = 
     General model:
     coeffsyDarkBg(x) = g+(1-g-l)*0.5*(1+erf((x-u)/sqrt(2*v^2)))
     Coefficients (with 95% confidence bounds):
       g =    0.009043  (-0.008105, 0.02619)
       l =     0.01047  (-0.05935, 0.08029)
       u =       7.492  (7.274, 7.71)
       v =        1.83  (1.584, 2.077)

coeffsyDarkBg.u % Bias
ans =
    7.4921

coeffsyDarkBg.v % Discrimination sensitivity

ans =

    1.8305

The units for these values are the same as for the x axis (see the previous post for details on their interpretation).

By way of comparison, the values obtained for the bias discrimination sensitivity for this data for the obscured target condtion using FitPsychcurveLogit were 5.5 and 2.1. Here they are

coeffsObscured = 

...

       u =         5.5  (5.325, 5.675)

       v =        1.94  (1.68, 2.2)

Even with such clean data, there's an noticeably better fit of the curve to the data, and a better estimate of the discrimination sensitivity parameter (which is a rather important parameter!) Obviously, we need to use statistics to really test the fit and find its reliability, but that's a topic for another time. For those interested, I'd suggest the two Wichmann and hill papers as further reading: 2001a, 2001b.

Code run-through
function [ffit, curve] = ...
    FitPsycheCurveWH(xAxis, yData, varargin)

FitPsychCurveWH requires two inputs the xaxis (xAxis) and the data to fit to (yData). In addition it can take a third input, which should contain the start points and upper and lower limits for the coefficients arranged in a 3x4 matrix where each column represents a coefficient (g, l, u, v, respectively). Row 1 contains the upper limits, row 2 the start points, and row 3 the lower limits. When this input is included it is caught by the function varargin, which stores the values in a variable also named varargin. The variable varargin is a cell array; if the matrix of limits is specified, it will be stored in varagin{1}, if the limits are not entered, varargin will be empty.

% Start points and limits
if numel(varargin{:})>1
    useLims=1;
    UL=varargin{1}(1,:);
    SP=varargin{1}(2,:);
    LM=varargin{1}(3,:);
else
    useLims=0;
end

Firstly the function checks the contents of varagin. If varagin isn't empty, it extracts the upper and lower limits and start points and stores them as vectors in their own variables. This is done mainly to make the coding below a bit clearer. The { } brackets extract the contents of the first cell of varargin, and then the ( ) brackets extract the relevant row and and columns from the matrix stored in varargin's cell. The variable useLimts will be used later to control the input to fit.

% Transpose if necessary
if size(xAxis,1)<size(xAxis,2)
    xAxis=xAxis';
end
if size(yData,1)<size(yData,2)
    yData=yData';
end

Next the input data and axis are transposed to columns if they've been entered as rows, to stop fit complaining.

% Check range of data
if min(yData)<0 || max(yData)>1  
     % Attempt to normalise data to range 0 to 1
    yData = yData/(mean([min(yData), max(yData)])*2);
end

The data is normalised to the range 0 to 1 if any values extend above 1 or below 0. This is done by finding the midpoint (ie. the mean of the max and min values) and then dividing the whole of yData by the midpoint * 2.
 
% Prepare fitting function
F=@(g,l,u,v,x) g+(1-g-l)*0.5*(1+erf((x-u)/sqrt(2*v^2)));

This is the definiation of Wichmann and Hill's function which is stored as an anonymous function in the variable F. "Anonymous" just refers to the fact the function is stored in a variable in the workspace, instead of in a .m file.

% Fit using fit function from fit toolbox
if useLims==1
    % Start points and limits specified, use while fitting
    ffit=fit(xAxis, yData, F, 'StartPoint', SP, 'Upper', UL, 'Lower', LM);
else
    % Fits not specified, don't use while fitting
    ffit=fit(xAxis, yData, F);
end

Now everything is prepared, the data and F are passed to Matlab's fit function. If useLims was set to 1 earlier, the start points and upper and lower limits are also passed to fit as parameters. fit outputs a fit object (called ffit here) containing the coefficients of the fit and their confidence bounds.

% Create a new xAxis with higher resolution
fineX = linspace(min(xAxis),max(xAxis),numel(xAxis)*50);
% Generate curve from fit
curve = feval(ffit, fineX);
curve = [fineX', curve];

Finally, a higher resolution curve is generated for plotting (and to provide an example of using feval!). feval is used to evaluate the fit object ffit over the points in the x axis specified (fineX). Alternatively, the fit object can also be passed directly to functions such as plot (eg. plot(ffit)).