PSYCHOMETRICA

Functions for Fitting and Plotting Psychometric Data

©1995-2018 Andrew B. Watson

Credits & Version

Copyright 1995-2018 by Andrew B. Watson. All rights reserved. Do not distribute or reproduce without permission. No warranties expressed or implied. This Version is compatible with Mathematica 10.

Date of this version:

DateString[]

History

1.0->2.0

Changed label location in PsychometricPlot.
Took shift argument out of all psychometric functions, and all examples using such functions.
Several arguments changed to options in PsychometricFit. Added FitLimits option.
PsychometricFit now accepts a list of tables or the output of QuestExperiment.
PsychometricPlot accepts multiple fits.
Added PNames to each psychometric function.
Complete print of parameters in PsychometricPlot.
Underflow messages turned Off and On in psychometric functions.
Term "value" changed to "strength" throughout.

2.0->2.1

added LinearPF 6/30/96

2.1->2.2

dropped LinearPF
added UniformPF
added ArgumentNames
dropped PNames data for each PF.
lightened gray of errorbars (.5-.75)
moved errorbars behind data points

2.2-3.0

Changed underflow message name from "under" to "unfl" in LogWeibullPF and NormalPF. 9/22/97

3.0->3.1

Improved automatic initial estimate of threshold in PsychometricFit. Option PlotError added.

Removed Statistics`ContinuousDistributions` from initializations.
Removed ArgumentNames (is in Miscellany.m).
Added Off[General::spell1]

Off[Power::"infy"] in PsychometricFit

{2004,9,13,16,39,24.535654}

Fixed PsychometricFit to work with new FindMinimum syntax.

{2007, 10, 16, 14, 58, 37.944963}

Modified LogWeibullPF to avoid overflow.

{2010, 2, 16, 17, 54, 29.131186}

Slight change to formatting of legend in PsychometricPlot,

{2012,7,19,9,50,31.664680}

Added lapse to CNormalPF

{2014,8,4,10,24,6.352587}

Many updates to plots, function definitions simplified to use CDF where possible.
Introduced WeibullPF.
Streamlined PsychometricFit
Streamline PsychometricPlot
Added FitMethod
Integrated inverse functions into main functions
Created PlotPF
Created TimingPF
Copied some functions from Miscellany and Quest (see utility functions)
Added trials to legend in PsychometricPlot

{2014,12,29,13,44,31.965248}

PsychometricFit: Changed calculation of error function. Added binomial error. Changed plot of likelihood so that it is now actual log likelihhod. Changed final reported error to agree with Chi-Square statistic as descibed in Watson (1979). Changed first FitLimit to Automatic, which means use range of actual strengths. Changed limits in optimization to params[[i]]+ {0,.001} to prevent error when params[[i]]==0.
LogWeibullPF: reverted to algebraic formula beacuse version using GumbelDistribution was blowing up under certain circumstances.
-ABW

{2015, 4, 18, 14, 52, 49.743290}

Made definition of MinMax conditional on its absence (i.e. in Mathematica versions 9 and below).

{2018, 3, 14, 11, 35, 57.637374}

jas Turned “General::munfl” off.

{2020,4,25,12,16,48.726277}

jas cleaned up PreciseErrorBars section.

Introduction

These functions implement general fitting and plotting of psychometric data. Weibull, Normal and Logistic psychometric functions are also defined. Some background reading on the fitting method is available in:

Watson, A. B. (1979). Probability summation over time. Vision Research 19, 515-522.

This Notebook begins with function definitions and some documentation. If you are not concerned with implementation, and want only to use the functions, take a look at the Tutorial.

Notes on implementation

In the current implementation we have tried where possible to rely on built-in probability distributions. This makes the implementation simpler and more transparent, but may also make it slower. However it is rare that speed is critical in calculation of psychometric functions. We have also used integers constants where possible, and have not coerced numerical results through use the the N function. This allows symbolic manipulation. Where speed is required, ensuring that parameters are real and  coercing to real may be helpful. As a last resort, a variant can be written with numerical expressions.

Preliminaries

PsychometricaStyles={Frame->True,BaseStyle->{12,FontFamily->"Helvetica"},ImageSize->{Automatic,240}};

Off[General::munfl]

Utility Functions

MinMax

If[Names["MinMax"]=={},MinMax::usage = "MinMax[list_]: Returns a list {min,max} of the largest and smallest entries in the list.";MinMax[list_List] := {Min[list],Max[list]}]

dB

From Miscellany.nb

The decibel (dB) is 1/20 of the log to the base 10. These functions convert to and from dB.

Clear[dB]

dB[x_] := 20 Log[10,x]

dBInverse

From Miscellany.nb

dBInverse[x_]:=((10 )^(x /20))

PlotPF

Clear[PlotPF]

PlotPF::usage = "PlotPF[function_,thlist_ ,slopelist_]: This utility function shows a plot of a specified psychometric function and allows the user to explore different parameters. The first argument is the name of the function. The second argument is a list {min,max,start) for threshold, and the third is a similar list for slope. The green arrows show the location of threshold, and the probability when strength equals threshold. The red arrows show the strength required to achieve a selected probability. This is obtained using the inverse psychometric function.";

PlotPF[fn_,thlist_ ,slopelist_] :=
 Manipulate[
(target0=fn[threshold,{threshold,slope,guess,lapse}];
strength0=fn[target,{threshold,slope,guess,lapse},True];
Plot[fn[x,{threshold,slope,guess,lapse}],{x,thlist[[1]],thlist[[2]]},PlotRange->{0,1}
,FrameLabel->{"Strength","Probability"}
,PlotLabel->ToString[fn]
,Evaluate@PsychometricaStyles
,Epilog->{Green
,Point[{threshold,target0}]
,Arrow[{{threshold,target0},{thlist[[1]],target0}}]
,Arrow[{{threshold,target0},{threshold,0}}]
,Text["threshold target = "<>ToString[NumberForm[N@target,2]],{thlist[[1]],target0},{-1.1,-2}]
,Text["threshold = "<>ToString[NumberForm[N@threshold,2]],{threshold ,.2},{-1.2,0}]
,Red
,Point[{strength0,target}]
,Arrow[{{strength0,target},{thlist[[2]],target}}]
,Arrow[{{strength0,target},{strength0,0}}]
,Text["target = "<>ToString[NumberForm[N@target,2]],{thlist[[2]],target},{1.1,-2}]
,Text["strength = "<>ToString[NumberForm[N@strength0,2]],{strength0 ,.1},{-1.2,0}]}])
,{{threshold,thlist[[3]]},thlist[[1]],thlist[[2]]}
,{{slope,slopelist[[3]]},slopelist[[1]],slopelist[[2]]}
,{{guess,0},0,.5}
,{{lapse,0},0,.5}
,{{target,target0},0,1}
,ContinuousAction->True]

PlotPF[LogisticPF,{0,1,.5},{1,10,4}]

TimingPF

Clear[TimingPF]

TimingPF[fn_,p_, min_,max_,n_:10000] := Block[{result},
{Timing[result=fn[N@Range[min,max,(max-min)/n],p]][[1]],
Total[result]}]

TimingPF[WeibullPF,{.5, 4, .1, .01},0,1]

IndexMin

Returns the index of the smallest element in a list.

From Miscellany.nb

IndexMin[list_] :=
 First[Flatten[Position[list, Min[list]]]]

IndexMax

IndexMax::usage = "IndexMax[list_]: This utility function returns the index of the largest element in a list (in the case of a tie, it returns the first of the winners).";

IndexMax[list_] := First[Flatten[Position[list,Max[list]]]]

ArgumentNames

Returns a list of lists of strings corresponding to the names of the arguments in the definitions of a function, one list per definition.

From Miscellany.nb

ArgumentNames[f_] :=( ToString /@ 
Function[patt,
  Cases[patt,
	  p_Pattern :> HoldForm @@ HeldPart[p, 1],
	   {0, Infinity},Heads->True]
	   ][#[[1]]])& /@ DownValues[f]

g[d_,{a_,b_,c_}] := 1

res1= ArgumentNames[g]

Psychometric Functions

The following are definitions for a number of commonly used psychometric functions, which describe probability of a positive response (yes, or correct) as a function of some measure of stimulus strength. They are typically cumulative probability distribution functions, ammended by guessing and lapsing factors. All the functions defined here follow a standard form f[x_,p_List] where p is a list of parameters. This standard form allows the functions to be used by fitting routines and by the Quest staircase procedure (Quest.nb). User-written functions should follow the same convention to obtain these capabilities. For many of the functions, the parameter list has the form {threshold, slope, guess, lapse}. Threshold is the strength at which some designated probability is obtained. Slope is a measure of the steepness of the function in the neighborhood of threshold. Guess, is the lower asymptotic probability, and lapse is the probability of an error at very high strength, that is, 1 minus the upper asymptotic probability.

WeibullPF

Introduction

The Weibull function has proven useful as a description of the human psychometric function for visual detection. The following are some useful references.

Nachmias, J. (1981). On the psychometric function for contrast detection. Vision Research 21(2), 215-223.

Quick, R. F. (1974). A vector magnitude model of contrast detection. Kybernetik 16, 65-67.

Watson, A. B. (1979). Probability summation over time. Vision Research 19, 515-522.

Weibull, W. (1951). A statistical distribution function of wide applicability. J. appl Mech. 18, 292-297.

Function

Clear[WeibullPF]

WeibullPF::usage = "WeibullPF[strength_, {threshold_, slope_, guess_, lapse_}, inverse_:False]: This is a psychometric function based on the Weibull probability distribution. The stimulus strength is x. Threshold is the strength at which some target probability is obtained. Slope is a measure of the steepness of the function in the neighborhood of threshold. Guess, is the lower asymptotic probability,and lapse is the probability of an error at very high strength, that is, 1 minus the upper asymptotics probability. If inverse is True, it computes the inverse psychometric function: given a probability p, it returns the stimulus strength.";

WeibullPF[strength_,{threshold_,slope_,guess_,lapse_},inverse_:False] :=  
 guess + (1-guess-lapse) CDF[WeibullDistribution[slope,threshold],strength]

WeibullPF[p_,{threshold_,slope_,guess_,lapse_},True] := threshold (((Log[((-1+guess+lapse)/(-1+lapse+p))])^(1/slope)))

Plot

PlotPF[WeibullPF,{0,1,.5} ,{1,10,4}]

Timing

TimingPF[WeibullPF,{.5, 4, .1, .01},0,1]

Derivation

CDF[WeibullDistribution[slope,threshold],x]

dBInverse@LogWeibullPFInverse[p,{dB[threshold],slope,guess,lapse}]//FullSimplify

LogWeibullPF

This is identical to the WeibullPF, except that strength and threshold are expressed in decibels (dB). It is equivalent to the Gumbel distribution. The LogWeibullPF function is also related to the Log Weibull distribution, the Extreme Value distribution, and the Fisher-Tippet distribution.

Clear[LogWeibullPF]

LogWeibullPF::usage = "LogWeibullPF[strength_,{threshold_,slope_,guess_,lapse_}]]: This is a psychometric function based on the Weibull probability distribution, where strength and threshold are defined in decibels.  Threshold is the strength at which some target probability is obtained. Slope is a measure of the steepness of the function in the neighborhood of threshold. Guess is the lower asymptotic probability,and lapse is the probability of an error at very high strength, that is, 1 minus the upper asymptotics probability. If inverse is True, it computes the inverse psychometric function: given a probability p, it returns the stimulus strength.";

LogWeibullPF[x_,{threshold_,slope_,guess_,lapse_},inverse_:False] :=  
guess + (1.-guess-lapse) (1.-Exp[- (dBInverse[x - threshold]) ^ slope])

LogWeibullPF[p_,{threshold_,slope_,guess_,lapse_},True] := threshold+((Log[((Log[((-1+guess+lapse)/(-1+lapse+p))])^(20/slope))])/(Log[10]))

Using GumbelDistribution

Unfortunately this blows up in timing.

Clear[LogWeibullPF1]

LogWeibullPF1[strength_,{threshold_,slope_,guess_,lapse_},inverse_:False] :=
guess+(1.-guess-lapse)CDF[GumbelDistribution[threshold,(20./(slope Log[10.]))],strength]

AbsoluteTiming[LogWeibullPF[14,{-40.,4.,0.5,.01}]]

AbsoluteTiming[LogWeibullPF1[14,{-40.,4.,0.5,.01}]]

Using ExtremeValueDistribution

The LogWeibullPF could be defined in terms of the ExtremeValueDistribution, like so:

LogWeibullPF2[x_,{threshold_,slope_,guess_,lapse_}] :=
guess+(1.-guess-lapse)(1.-CDF[ExtremeValueDistribution[threshold,20./(slope Log[10.])],-x])

Derivation

Clear[slope]

CDF[WeibullDistribution[slope,threshold],x]

CDF[ExtremeValueDistribution[a,b],y]

Exp[-dBInverse[dB[x] - dB[threshold]]^slope]//Simplify

E^(-(E^((-Log[threshold] + Log[x])))^slope)

E^(-(E^(slope (-Log[threshold] + Log[x]))))

CDF[WeibullDistribution[slope,dBInverse[threshold]],dBInverse[x]]

((((10^(-(threshold/20)+(x/20)))))^slope)//Simplify

CDF[GumbelDistribution[a,b],y]

10^Log[10,E]//Simplify

Solve[-((slope(threshold-x))/20)==Log[10,E]((-a+y)/b),{a,b,y}]

Solve[slope/20 == Log[10,E]/b,b]

Solve[slope/20 == Log[10,E]/b,b]

CDF[GumbelDistribution[threshold,(20/(slope Log[10]))],x]

((10^(slope((-threshold+x)/20))))==(10^((1/20) slope (-threshold+x)))

Using WeibullDistribution

LogWeibullPF3[x_,{threshold_,slope_,guess_,lapse_}] :=  
 guess + (1-guess-lapse) CDF[WeibullDistribution[slope,dBInverse[threshold]],dBInverse[x]]

Compare alternate versions

This returns the timing and summed values.

Clear[try]

try[fn_,limit_:60] := 
{First@AbsoluteTiming[tst=Table[fn[y, {0., 4., 0., 0.}],{y,-limit,limit,.1}];],
Plus @@ Flatten[tst]//N}

try[LogWeibullPF]

try[LogWeibullPF1,30]

try[LogWeibullPF1,40]

This takes too long:

try[LogWeibullPF1,50]

try[LogWeibullPF2]

try[LogWeibullPF3]

try /@ {LogWeibullPF,LogWeibullPF2,LogWeibullPF3}//TableForm

Bounds (not used)

Considering possible bounds on computation. Beyond about {-60,20}/beta the results are constant.

LogWeibullPF[20/b,{0,b,0,0}] //PowerExpand//N

LogWeibullPF[-60/b,{0,b,0,0}] //PowerExpand//N

beta=2;lims={-60,20}/beta;
LogWeibullPF[#,{0,beta,0,0}] & /@ lims//N
Plot[LogWeibullPF[x,{0,beta,0,0}],{x,-30,30}
,Epilog->{Line[{{ #,0},{ #,1}}]& /@ lims}]

Example

LogWeibullPF[Range[-6,6],{0.,4.,0,.01}]

Plots

PlotPF[LogWeibullPF,{-6,6,0} ,{1,10,4}]

beta = 2^Range[0,4];
functions=LogWeibullPF[x,{0.,#,0,0.01}]& /@ beta;
Plot[functions,{x,-6,6},PlotRange->{0,1}
,FrameLabel->{"Strength","Probability"}
,PlotLegends->Placed[LineLegend[beta,LegendLabel->"beta"],{Right,Bottom}]
,Evaluate@PsychometricaStyles]

guess ={0.5`,0.25`,0};
functions=LogWeibullPF[x,{0.,4,#,0.01}]& /@ guess;
Plot[functions,{x,-6,6},PlotRange->{0,1}
,FrameLabel->{"Strength","Probability"}
,PlotLegends->Placed[LineLegend[guess,LegendLabel->"guess"],{Right,Bottom}]
,Evaluate@PsychometricaStyles]

lapse ={0,.1,.2};
functions=LogWeibullPF[x,{0.,4,0,#}]& /@ lapse;
Plot[functions,{x,-6,6},PlotRange->{0,1}
,FrameLabel->{"Strength","Probability"}
,PlotLegends->Placed[LineLegend[lapse,LegendLabel->"lapse"],{Right,Bottom}]
,Evaluate@PsychometricaStyles]

Comparison with Weibull

The Weibull PF is identical to the LogWeibull exept that the latter has logarithmic (dB) values for x and threshold.

LogWeibullPF2[x_,{threshold_,slope_,guess_,lapse_}] :=  
 guess + (1-guess-lapse) CDF[WeibullDistribution[slope,dBInverse[threshold]],dBInverse[x]]

CDF[WeibullDistribution[slope,threshold],x]

{threshold,slope,guess,lapse}= {0,4,.5,.01};
Plot[{LogWeibullPF[x,{threshold,slope,guess,lapse}]
,LogWeibullPF2[x,{threshold,slope,guess,lapse}]
,WeibullPF[dBInverse[x],{dBInverse[threshold],slope,guess,lapse}]}
,{x,-6,6},PlotRange->{0,1}
,FrameLabel->{"Strength","Probability"}
,PlotStyle->{Red,Dashed,Dotted}
,Evaluate@PsychometricaStyles]

Derivation of inverse

Solve[guess+(1-guess-lapse) (1-Exp[-((dBInverse[x-threshold])^slope)])==p,x]

Relation to PsychToolbox

By default PsychToolbox’s QUEST routine uses PelliPF[x,{threshold,beta,gamma,delta}]=LogWeibullPF[x,{threshold,20*beta,gamma,delta*(1-gamma)}]

LogisticPF

Clear[LogisticPF]

LogisticPF::usage = "LogisticPF[strength_, {threshold_, slope_, guess_, lapse_}]: This is a psychometric function based on the Logistic probability distribution.  Threshold is the strength at which some target probability is obtained. Slope is a measure of the steepness of the function in the neighborhood of threshold. Guess is the lower asymptotic probability,and lapse is the probability of an error at very high strength, that is, 1 minus the upper asymptotics probability. If inverse is True, it computes the inverse psychometric function: given a probability p, it returns the stimulus strength.";

LogisticPF[strength_,{threshold_,slope_,guess_,lapse_},inverse_:False]:=  
 guess + (1 - guess - lapse) CDF[LogisticDistribution[threshold,1/slope],strength]

LogisticPF[p_,{threshold_,slope_,guess_,lapse_},True] := ((slope threshold+Log[((guess-p)/(-1+lapse+p))])/slope)

Example

LogisticPF[Range[0,1,.1],{.5,4.,0.1,.01}]

Derivation

CDF[LogisticDistribution[μ,β],x]

Plot

PlotPF[LogisticPF,{0,1,.5},{1,10,4}]

Timing

Retina Pro 2.6 GHz Intel Core i7

TimingPF[LogisticPF,{.5,4,.5,.01},0,1]

Old

LogisticPF[x_,{threshold_,slope_,guess_,lapse_}]:=  
 N[guess + (1 - guess - lapse) /
   (1 + Exp[ slope (threshold - x)])]

Derivation of invserse

Solve[{LogisticPF[x,{threshold,slope,guess,lapse}]==p,x>0},x,Reals]

NormalPF

Clear[NormalPF]

NormalPF::usage = "NormalPF[strength_,{threshold_,slope_,guess_,lapse_}]: The NormalPF psychometric function is the cumulative distribution function of a Normal density with mean threshold and standard deviation 1/slope. The stimulus strength is x. The guessing rate is guess, the lapsing rate is lapse. If inverse is True, it computes the inverse psychometric function: given a probability p, it returns the stimulus strength.";

NormalPF[strength_,{threshold_,slope_,guess_,lapse_},inverse_:False] := guess+(1-guess-lapse) CDF[NormalDistribution[threshold,1/slope],strength]

NormalPF[p_,{threshold_,slope_,guess_,lapse_},True] := threshold-(Sqrt[2]) (1/slope) InverseErfc[((2 (guess-p))/(-1+guess+lapse))]

Example

NormalPF[Range[0,1,.1],{.5,4.,0.1,.01}]

Plots

PlotPF[NormalPF,{0,1,.5},{1,10,4}]

Timing

Retina Pro 2.6 GHz Intel Core i7

TimingPF[NormalPF,{.5,4,.1,.01},0,1]

Derivation

CDF[NormalDistribution[threshold,1/slope],strength]

Old

NormalPF[x_,{threshold_, slope_, guess_, lapse_}] :=   
 Module[{p}, 
  Off[General::unfl];
  p = guess + 0.5 (1. - guess - lapse) *
   (1.+ Erf[(0.7071067811865476 slope (x-threshold))]);
 On[General::unfl];
 p]

NormalPF[x_,{threshold_, slope_, guess_, lapse_}] :=   
guess + 0.5 (1. - guess - lapse) *
   (1.+ Erf[(0.7071067811865476 slope (x-threshold))])

Derivation of inverse

Solve[NormalPF[x,{threshold,slope,guess,lapse}]==p,x]

CNormalPF

This is identical to the NormalPF, except that uses sd instead of slope.

Clear[CNormalPF]

CNormalPF::usage = "CNormalPF[strength_, {mean_, sd_, guess_, lapse_}]: The CNormalPF psychometric function is the cumulative distribution function of a Normal density with mean mean and standard deviation sd. The stimulus strength is strength. The guessing rate is guess, the lapsing rate is lapse. The function returns a probability of success.";

CNormalPF[strength_,{mean_,sd_,guess_,lapse_},inverse_:False] := NormalPF[strength,{mean,1/sd,guess,lapse}]

CNormalPF[p_,{mean_,sd_,guess_,lapse_},True] := mean-(Sqrt[2]) sd InverseErfc[((2 (guess-p))/(-1+guess+lapse))]

CNormalPFInverse[p_,{mean_,sd_,guess_,lapse_}] := CNormalPF[p,{mean,sd,guess,lapse},True]

UniformPF

Clear[UniformPF]

UniformPF::usage = "UniformPF[strength_,{threshold_,slope_,guess_,lapse_}]: The UniformPF psychometric function is based on the cumulative uniform distribution of width 1/slope, centered at threshold, ammended by guessing and lapsing parameters. It represents a linear increase in probability from a lower asymptote of quess to an upper asymptote of 1-lapse, with a rate of increase of (1 - guess - lapse)*slope. At threshold, the probability is (1 + guess - lapse)/2 .  If inverse is True, it computes the inverse psychometric function: given a probability p, it returns the stimulus strength.";

UniformPF[strength_,{threshold_, slope_, guess_, lapse_},inverse_:False] :=
guess + (1 - guess - lapse) CDF[UniformDistribution[threshold +{-1,1}/2/slope],strength]

UniformPF[p_,{threshold_, slope_, guess_, lapse_},True] :=
((1+guess-lapse-2 p)/(2 (-1+guess+lapse) slope))+threshold

Plot

PlotPF[UniformPF,{0,1,.5},{1,10,4}]

Derivation of inverse

Solve[UniformPF[x,{threshold,slope,guess,lapse}]==p,x]

UniformPF[x,{threshold,slope,guess,lapse}]

Solve[guess + (1 - guess - lapse)(slope ((1/(2 slope))+strength-threshold))==p,strength]

((1+guess-lapse-2 p-2 slope threshold+2 guess slope threshold+2 lapse slope threshold)/(2 (-1+guess+lapse) slope))//FullSimplify

Old

UniformPF[strength_,{threshold_, slope_, guess_, lapse_}] :=
guess + (1. - guess - lapse) *
 (Max[0,Min[1,slope (strength-threshold+1/(2 slope))]])

PsychometricFunction

Not yet modified to include invser functions.

Clear[PsychometricFunction]

PsychometricFunction::usage = "PsychometricFunction[strength_, {threshold_, slope_, guess_, lapse_}, opts___Rule]: The psychometric function defines the probability of being correct in a psychophysical task as a function of signal strength. The function is based on the cumulative probability distribution, modified by guessing and lapsing probabilities. The distribution can be specified by the Distribution option, which degaults to WeibullDistribution.   Threshold is the strength at which some target probability is obtained. Slope is a measure of the steepness of the function in the neighborhood of threshold. Guess, is the lower asymptotic probability, and lapse is the probability of an error at very high strength, that is, 1 minus the upper asymptotics probability. The option StrengthFunction specifies how to transform from strength and threshold specified in some new units to the standard units. For example, if strength and threshold are specified in dB, the use StrengthFunction->dBInverse.";

Options[PsychometricFunction] = {Distribution->WeibullDistribution,StrengthFunction->Identity};

PsychometricFunction[strength_,{threshold_, slope_, guess_, lapse_},opts___Rule] := Block[{distribution,parameters,sfunction},
{distribution,sfunction} = {Distribution,StrengthFunction} /. {opts} /. Options[PsychometricFunction];
parameters =
Switch[distribution
,WeibullDistribution,Sequence@@{slope,sfunction@threshold}
,NormalDistribution,Sequence@@{sfunction@threshold,1/slope}
,LogisticDistribution,Sequence@@{sfunction@threshold,1/slope}
,UniformDistribution,sfunction[threshold +{-1,1}/2/slope]];
guess + (1 - guess - lapse) CDF[distribution[parameters],sfunction[strength]]
]

Examples

NormalPF[.6,{.5, 4, .1, .01}]

PsychometricFunction[.6,{.5, 4, .1, .01},Distribution->NormalDistribution]

PsychometricFunction[dB[.6],{dB[.5], 4, .1, .01},Distribution->NormalDistribution,StrengthFunction->dBInverse]

Timing

TimingPF[WeibullPF,{.5, 4, .1, .01},0,1]

TimingPF[PsychometricFunction,{.5, 4, .1, .01},0,1]

Fitting and Plotting Psychometric Data

PsychometricFit

PsychometricFit::usage = "PsychometricFit[data_]: Fits a set of psychometric data with a psychometric function. The function returns the input data, the function name, and the specified and estimated parameters, and the residual error. The format of the returned material is {data, function, {parameters, error}}. The residual error is -2 times the log of the ratio of constrained and unconstrained hypotheses, which is asymptotically Chi-Square. The fit is based upon a maximum likelihood method, as described in: Watson, A. B. (1979) Probability summation over time, Vision Research 19, 515-522. data is a list of triples of the form {strength,no,yes}, where strength is a measure of signal strength, no is the number of no or incorrect responses, and yes is the number of yes or correct reponses. The function also accepts a list of such tables, and produces one fit per table. The function also accepts the data structure produced by QuestExperiment from the notebook Quest.nb. 
Options:
PsychometricFunction: name of the psychometric function to be fit (default = LogWeibullPF). The function must adhere to a standard format: function[x_, parameters_List]. 
PsychometricParameters: initial or final values of psychometric function parameters (default = {Automatic, 4., 0.5, .01}). The value Automatic in the first position indicates that the initial value will be estimated from the data. This estimate is obtained by evaluating the log likelihood function, using the default parameters for the psychometric function, at each intensity used, and selecting the intensity with the largest likelihood. This log likelihood function can be plotted by selecting the option PlotLikelihood->True. 
FitParameters: specifies whether each parameter should be estimated (default = {True}). Unspecified values default to False. For example, to estimate only the first and third parameters, use FitParameters->{True,False,True}. By convention, the first parameter is threshold and the second parameter is slope. This convention is relevant only in functions such as PsychometricPlot, which accepts directly the output of PsychometricFit. 
FitLimits: bounds on each parameter (default = {Automatic, {0,$MaxMachineNumber}, {0,1}, {0,1}}
FitMethod: the function used to optimize parameters (default = FindMinimum). As an alterantive the method NMinimize may be specified, which returns a more global optimum.";

Options[PsychometricFit] = {PsychometricFunction->LogWeibullPF, PsychometricParameters->{Automatic,4.,0.5,0.01}, FitParameters->{True},
FitLimits->{Automatic,{0,$MaxMachineNumber},{0,1},{0,1}},
PlotLikelihood->False,FitMethod->FindMinimum,Verbose->False};

PsychometricFit[data_,opts___Rule]:=Block[{fit,function,limits,params,err,args,strength,no,yes,total,result,dummy,dum1,dum2,dum3,dum4,dum5,dum6,plotlike,fmethod,m0,m1,df,logbinomial,terr,terrf,verbose},{fit,params,function,limits,plotlike,fmethod,verbose}={FitParameters,PsychometricParameters,PsychometricFunction,FitLimits,PlotLikelihood,FitMethod,Verbose}/. {opts}/. Options[PsychometricFit];
{strength,no,yes}=Transpose[data];total=no+yes;dummy=Take[{dum1,dum2,dum3,dum4,dum5,dum6},Length[params]];

fit=Join[fit,Table[False,{Length[params]-Length[fit]}]];
df = Length[total]-Count[fit,True];
m1:=Total[NLogP[yes,#1]+NLogP[no,(1-#1)]]&[function[strength,dummy]];
logbinomial =N@Total[Log[Binomial[total,yes]]];
m0=N@Total[NLogP[yes,yes/total]+NLogP[no,no/total]];

terr=m1/. MapThread[#1->#2&,Rest/@{dummy,params}];
terrf=logbinomial+(terr/.{dum1->#1}&)/@strength;
If[plotlike,Print @ListLinePlot[Transpose[{strength,terrf}]
,FrameLabel->{"Strength","Log Likelihood"},Mesh->True,PsychometricaStyles]];
If[params[[1]]===Automatic,params[[1]]=strength[[IndexMax[terrf]]]];
If[limits[[1]]=== Automatic,limits[[1]]=MinMax[strength]];
If[verbose,Print["limits =",limits,"\ndf = ",df,"\nlogbinomial = ",logbinomial]];

result=Switch[fmethod,
FindMinimum,
(args={};
Table[If[fit[[i]]
,args=Append[args,Flatten[{dummy[[i]],params[[i]]+ {0,.001},limits[[i]]}]]
,dummy[[i]]=params[[i]]],{i,Length[params]}];
args=Prepend[args,-m1];
Check[FindMinimum@@args,{Null,(#1->Null&)/@dummy},{FindMinimum::"fmlim",FindMinimum::cvmit}]),
NMinimize,
(vars={};
cons = {};
Table[If[fit[[i]]
,(vars=Append[vars,dummy[[i]]];
cons = Append[cons,limits[[i,1]]<dummy[[i]]<limits[[i,2]]])
,dummy[[i]]=params[[i]]],{i,Length[params]}];
args={Prepend[cons,-m1],vars};

Check[NMinimize@@args,{Null,(#1->Null&)/@dummy}])];
{data,function,{dummy/. result[[2]],2(m0+result[[1]])}}]

PsychometricFit[data_,opts___Rule] :=
PsychometricFit[SequenceToTable[First[data]],opts] /;
 Depth[data] == 4 &&
 Not[Equal @@ (Last /@ (Dimensions /@ data))]

PsychometricFit[data_,opts___Rule] :=
 PsychometricFit[#,opts]& /@ data /;
 Depth[data] == 4 &&
 Equal @@ (Last /@ (Dimensions /@ data))

PsychometricPlot

PsychometricPlot::usage = "PsychometricPlot[fitdata]: Plots a set of psychometric data and a psychometric function, as returned by PsychometricFit. See that function for details. Also draws approximate Bernoulli confidence intervals around each data point, unless the option PsychErrorBars->False is given. To suppress printing of parameters and error of the fit, use the option PsychLegend->False. Options to ListPlot may also be included.";

Clear[PsychometricPlot]

Options[PsychometricPlot] = {PsychErrorBars->True, PsychLegend->True};

Formatted legend better abw 2/16/2010

added PlotRange->All to pfunc. 10/23/13.

PsychometricPlot[{data_,function_,{params_,error_}},opts___Rule]:=Module[{val,no,yes,vrange,vminmax,xminmax,pdata,pfunc,pbars,plegend,bars,legend,trials},{bars,legend}={PsychErrorBars,PsychLegend}/. {opts}/. Options[PsychometricPlot];{val,no,yes}=Transpose[data];vminmax={Min@#,Max@#}&@val;
vrange=-Subtract@@vminmax;
xminmax=vminmax+.1vrange{-1,1};
trials = Total[no+yes];
pbars=Graphics[{GrayLevel[0.75],MapThread[Line[Transpose[{{#1,#1},BinomialErrorBars[#2,#3]}]]&,{val,no,yes}]}];

plegend=Graphics[Text[Style[#,TextAlignment->Left]& @ 
(StringJoin@@Flatten[{ToString[function],Table[{"\n",ArgumentNames[function][[1,i+1]]," = ",ToString@NumberForm[params[[i]],3]},{i,Length[params]}],"\nerror = ",ToString@NumberForm[error,3]
,"\ntrials = ",ToString@trials}]),Scaled[{0,1}],{-1.2,1.1}]];

pdata=ListPlot[Transpose[{val,(yes/(yes+no))}]
,Sequence@@FilterRules[Flatten[{opts}],Options[ListPlot]]
,PlotRange->{xminmax,{-.02,1.02}}
,FrameLabel->{"Strength","Probability"}
,Axes->False
,PlotStyle->{PointSize[0.03]}
,PsychometricaStyles];
pfunc=Plot[function[x,params],{x,xminmax[[1]],xminmax[[2]]}
,PlotRange->All
,Evaluate[Sequence@@FilterRules[Flatten[{opts}],Options[Plot]]]
,PlotStyle->AbsoluteThickness[1]
,Evaluate @PsychometricaStyles];
Show[pdata,pfunc,If[legend,plegend,{}],If[bars,pbars,{}],pdata]]

fit1 = {{{2,4,3},{4,3,4},{6,10,34},{8,1,18},{10,0,10}},LogWeibullPF,{{6.373449280907232`,4.`,0.5`,0.01`},0.48960216028935327`}};

PsychometricPlot [fit1]

PsychometricPlot[data_,opts___Rule] :=
 PsychometricPlot[#,opts]& /@ data /;
  Depth[data] == 5

Examples

fit1 = {{{2,4,3},{4,3,4},{6,10,34},{8,1,18},{10,0,10}},LogWeibullPF,{{6.373449280907232`,4.`,0.5`,0.01`},0.48960216028935327`}};

PsychometricPlot[fit1]

PsychometricPlot[fit1
, PlotStyle -> {Black,PointSize[.04],AbsoluteThickness[3]}]

SequenceToTable

We use Block rather than Module because the latter does not seem to isolate pattern names.

SequenceToTable::usage = "SequenceToTable[QSequence]: Converts a QSequence data structure as produced by Quest into a table of the form {{{strength, #wrong, #right},...},...} where there is one list of triples for each condition. This is the data structure expected by PsychometricFit.";

SequenceToTable[seq_] := Block[{a,x,b,c,d,e,y},
 Map[Rest,
(Sort[#/.{{a__,0} ->{a,1,0},{a__,1} ->{a,0,1}}] //. 
{x___,{a__,b_,c_},{a__,d_,e_},y___} -> {x,{a,b+d,c+e},y})& /@
 Table[Cases[seq,{cond,__}]
 ,{cond,Max[First /@ seq]}] ,{2}]]

Approximate error bars

BinomialErrorBars::usage = "BinomialErrorBars[failures,sucesses]:  Compute approximate confidence intervals on the probability p of a  binomial random variable, given numbers of failures and successes.";

These functions are used to compute approximate confidence intervals for a binomial random variable (from Abramowitz & Stegun). An alternative solution using the statistical packages is described in the next section.

rsign[x_] := Floor[x / (x + 1.) + .6]

tail[no_,yes_] := Exp[Log[.05] / (yes + no)]

aprox1[no_,yes_] := (yes + .5 + 1.92 + 1.96 Sqrt[(yes + .5)*
 (1. - (yes + .5) /(yes + no)) + .96]) / (yes + no + 3.842)

aprox2[no_,yes_] := (yes - .5 + 1.92 - 1.96 Sqrt[(yes - .5)*
 (1. - (yes - .5) / (yes + no)) + .96]) / (yes + no + 3.842)

errbr1[no_,yes_] := rsign[no] * aprox1[no,yes] * rsign[yes] +
 (1. - rsign[no]) + (1. - rsign[yes]) * (1. - tail[no,yes])

errbr2[no_,yes_] := rsign[no] * aprox2[no,yes] * rsign[yes] +
 (1. - rsign[no]) * tail[no,yes]

BinomialErrorBars[no_,yes_] := {errbr2[no,yes],errbr1[no,yes]}

BinomialErrorBars[3,4]

Precise error bars

This section describes the construction of Binomial confidence limits by way of probability distributions. It can be slow, and it depends upon a root-finding operation which sometimes fails to converge.

We first write a likelihood function that describes the likelihood of obtaining a particular set of data {no,yes}, given a Binomial random variable with probability p.

lik0[no_,yes_,p_] := PDF[BinomialDistribution[yes+no,p],yes]

We expand this function.

lik0[no,yes,p]

We remove the conditionals (which are always true), defining a new function.

lik1[no_,yes_,p_] =
lik0[no,yes,p][[1,1,1]]

To determine the normalized pdf, we integrate lik1 to find the normalizing factor.

Integrate[ lik1[no,yes,p] ,{p,0,1}]

norm[no_,yes_] :=
 (Binomial[no + yes, yes] Gamma[1 + no] Gamma[1 + yes]) /
 Gamma[2 + no + yes]

We create a posterior pdf.

pdf[no_,yes_,p_] := lik1[no,yes,p] / norm[no,yes]

We will need the cdf of this density, so we get the indefinite integral.

Integrate[pdf[no,yes,p],{p,0,x}]//InputForm

Using this result, we define a cumulative distribution.

cdf[no_,yes_,x_] := (Beta[x, 1 + yes, 1 + no]*Gamma[2 + no + yes])/
  (Gamma[1 + no]*Gamma[1 + yes])

The confidence limit on the binomial probability p is reached when a probability alpha is reached in the cdf.

Clear[BinomialCL]

BinomialCL[no_,yes_,alpha_,opts___] := FindRoot[
   (Beta[x, 1 + yes, 1 + no]*Gamma[2 + no + yes])/
  (Gamma[1 + no]*Gamma[1 + yes]) == alpha
      ,{x,N[yes/(yes+no)]}
      ,opts][[1,2]]

This plot compares the approximate error bars defined in the previous section with the more precise ones defined here. The parameter n is the possible number of yes or no responses. Red points are upper limits, green are lower limits. The color saturation corresponds to percent correct.

n=20;Show[Graphics[Table[MapThread[{Hue[#1,(yes/(no+yes)),1],Point[{BinomialCL[no,yes,#2,AccuracyGoal->4,PrecisionGoal->4],#3[no,yes]}]}&,{{1,0.29`},{0.975`,0.025`},{errbr1,errbr2}}],{no,1,n},{yes,1,n}]],Frame->True,FrameLabel->{"Precise","Approximate"},PlotRange->{{0,1},{0,1}},Axes->False,AspectRatio->1]

Timing[res = Apply[Function[{x,y},BinomialCL[x,y,#,AccuracyGoal->4,PrecisionGoal->4]& /@ {.025,.975}],tst,1];]

NLogP

In the computation of log likelihoods, it is occasionally necessary to make use of the following rules, so that the log function does not blow up.

NLogP[0,0] := 0

NLogP[n_,0] := -Sign[n] Infinity

NLogP[n_,p_] := n Log[p]

SetAttributes[NLogP,Listable]

NLogP[0_,Indeterminate] := 0

NLogP[0,0.] := 0

Tutorial

In this section I provide some simple examples of how to use the functions PsychometricFit and PsychometricPlot, and how to create and use a new psychometric function.

Psychometric functions

A psychometric function is a theoretical relationship between the strength of a sensory signal (e.g. the contrast of a visual pattern) and the probability of success in a sensory task related to that signal (e.g. detecting the pattern). In general, as the signal gets stronger, the probability of success increases.

Several mathematical functions have been found to provide good descriptions of human psychometric data. In general, these functions tend to be cumulative probability distributions, modified by so-called guessing and lapsing parameters. The complete function is often parameterized by four numbers: the threshold, the slop, the guessing probability, and the lapsing probability.  Threshold is the strength at which some target probability is obtained. Slope is a measure of the steepness of the function in the neighborhood of threshold. Guess is the porbability of success at zero or very low strength (the lower asymptotic probability), and lapse is the probability of an error at very high strength (1 minus the upper asymptotics probability).

Consider the NormalPF:

?NormalPF

We can compute one value :

NormalPF[.6,{.5,2,0,0}]

We can compute several values :

NormalPF[Range[.1,1,.1],{.5,2,0,0}]

We can plot it and explore its parameters using the utility function PlotPF.

?PlotPF

PlotPF[NormalPF,{0,1,.5},{.1,10,4}]

Here is a a different psychometric function.

PlotPF[WeibullPF,{0,1,.5},{.1,10,1}]

Inverse Psychometric functions

it is often useful to have an inverse psychometric function, that given a probability yields a strength. This is useful, for example, if one would like to know the the strength required for some criterion probability (e.g. 0.75) but the target probability of the function is some other value.

We have provided inverse functions for several of the psychometric functions, called by setting the thrid argument to True. Here is an example:

WeibullPF[.75,{1,1,.5,.01},True]

Creating  a new psychometric function

Here we show how the user can create a new psychometric function. A psychometric function is easily created from any probabilty distribution. As an example, we choose the Cauchy distribution, which has infinite variance, and thus reminds us of certain observers we have known.

For two parameter distributions, the parameters are usually a location measure and a dispersion measure. The former is associated with the threshold, and the latter with the inverse slope. Sometimes the slope is equal to the dispersion measure (e.g. WeibullDistribution). Make sure you invert, or not, as required.

?CauchyDistribution

CauchyPF[strength_,{threshold_,slope_,guess_,lapse_},invserse_:False] :=  guess + (1 - guess - lapse) CDF[CauchyDistribution[threshold,1/slope],strength]

We plot and explore the paramaters.

Plot[CauchyPF[strength,{1,4,0,0}],{strength,0,2},Evaluate@PsychometricaStyles]

To create the inverse function, we attempt to solve for strength

Solve[CauchyPF[strength,{threshold,slope,guess,lapse}]==p,strength]

Then we write an appropriate function.

CauchyPF[p_,{threshold_,slope_,guess_,lapse_},True] :=  ((slope threshold+Tan[((2.45850922`*^8+2.45850922`*^8 guess-2.45850922`*^8 lapse-4.91701844`*^8 p)/(-1.56513558`*^8+1.56513558`*^8 guess+1.56513558`*^8 lapse))])/slope)

Now we can use the PlotPF function (which requires an inverse).

PlotPF[CauchyPF,{0,1,.5},{1,100,4}]

Fitting a psychometric function to a set of data

First we create some fictional data. These are triples of the form {strength, # wrong, # right}.

data1={{2,4,3},{4,3,4},{6,10,34},{8,1,18},{10,0,10}};

We fit the data. The function returns the data, the name of the fitted psychometric function, and a list containing a list of function parameters and a measure of the error of the fit.

fit1 = PsychometricFit[data1]

This result can be plotted directly using the function PsychometricFit, about which we will say more in a moment.

PsychometricPlot[fit1]

By default, only the first parameter (threshold) is estimated. To estimate both the threshold and slope, we use the option FitParameters. Note the reduction in the error.

fit2=PsychometricFit[data1,FitParameters->{True,True}]
PsychometricPlot[fit2]

The default LogWeibullPF psychometric function has four parameters (threshold, slope, guess, and lapse). Here we estimate threshold and guess, and we alter the default value for the slope parameter.

PsychometricPlot[PsychometricFit[data1,FitParameters->{True,False,True}
,PsychometricParameters->{Automatic,3.,0.5,0.01}]]

Fitting a different psychometric function

Another widely used psychometric function is the LogisticPF. To change the fitted function, we use the option PsychometricFunction.

PsychometricPlot[PsychometricFit[data1 ,PsychometricFunction->LogisticPF
  ,FitParameters->{True,True} ]]

Avaliable pschometric functions include LogWeibullPF, WeibullPF, LogisticPF, NormalPF, UniformPF.

Fitting and plotting several sets of data

Both PsychometricFit and PsychometricPlot accept multiple data sets, as when an experiment contains several conditions.  In this example each of two conditions contributes one list of triples.

data2 = {
  {{3, 3, 4}, {5, 10, 34}, {7, 1, 18} }
 ,{{3, 4, 3}, {5, 16, 34}, {7,5,37}, {9, 1, 18} } };

fit2 =PsychometricFit[data2]

PsychometricPlot[fit2]

What happens when a fit cannot be found?

If a fit cannot be produced, a structure of the correct shape is returned, but with Nulls in place of all parameters and the error (beeps and messages will also be produced). This sometimes happens when the parameter limits are exceeded, or when the error function has no minimum.

PsychometricFit[{{2,4,3}}]

Some plotting options

The basic plot includes proportion correct (points), fitted function (curve), binomial confidence limits (grey vertical lines), and printed parameters. The latter two are optional.

PsychometricPlot[fit1,PsychErrorBars->False]

PsychometricPlot[fit1,PsychLegend->False]

In addition, options for ListPlot may be included.

PsychometricPlot[fit1
,BaseStyle->{14,AbsoluteThickness[4],FontFamily->"Helvetica",FontWeight->Bold,FontSlant->Italic,FontColor->Blue}
,PlotStyle->{Red,Dashed,AbsoluteThickness[2],AbsolutePointSize[20]},Background->Yellow
,FrameLabel->None]

Problems and solutions

Local mimima

Using the default FitMethod->FindMinimum, the function PsychometricPlot may not return a global minimum.

PsychometricPlot[PsychometricFit[data1  ,PsychometricFunction->NormalPF
  ,FitParameters->{True,True}]]

This can be corrected by a better guess for the slope.

PsychometricPlot[PsychometricFit[data1
  ,PsychometricFunction->NormalPF
  ,FitParameters->{True,True}
,PsychometricParameters->{Automatic,1.,0.5,0.01}]]

It can also be corrected by using FitMethod->NMinimize.

PsychometricPlot[PsychometricFit[data1
  ,PsychometricFunction->NormalPF
  ,FitParameters->{True,True}
,FitMethod->NMinimize
 ]]

Development

These are functions under development.

GeneralPF

Not ready for prime time. The idea was to allow specification of an arbitrary two-parameter distribution in terms of mean and variance. The problem is that the “natural” parameters must be solved for. Makes it slow.

GeneralPF[distribution_,{mean_,variance_,guess_,lapse_}] := guess+(1-guess-lapse) CDF[distribution[Sequence @@ Drop[parameters,-2]],x]

GeneralPF[distribution_,{mean_,variance_,guess_,lapse_}] := Block[{},
dmean = Mean[distribution[a,b]];
dvar = Variance[distribution[a,b]];

guess+(1-guess-lapse) CDF[distribution[Sequence @@ Drop[parameters,-2]],x]

Plot[GeneralPF[NormalDistribution,{0,1,.1,.05}],{x,-3,3},PlotRange->{0,1}
,FrameLabel->{"Strength","Probability"}
,Evaluate@PsychometricaStyles]

Mean[WeibullDistribution[a,b]]

Variance[WeibullDistribution[a,b]]

Clear[mean]

Solve[{mean == Mean[WeibullDistribution[a,b]],variance == Variance[WeibullDistribution[a,b]],mean>0,variance>0,a>0,b>0},{a,b}]

mean=0.5;
variance=.5;
NSolve[{mean == Mean[WeibullDistribution[a,b]],variance == Variance[WeibullDistribution[a,b]],a>0,b>0},{a,b},Reals]