**QUEST+ is a Bayesian adaptive psychometric testing method that allows an arbitrary number of stimulus dimensions, psychometric function parameters, and trial outcomes. It is a generalization and extension of the original QUEST procedure and incorporates many subsequent developments in the area of parametric adaptive testing. With a single procedure, it is possible to implement a wide variety of experimental designs, including conventional threshold measurement; measurement of psychometric function parameters, such as slope and lapse; estimation of the contrast sensitivity function; measurement of increment threshold functions; measurement of noise-masking functions; Thurstone scale estimation using pair comparisons; and categorical ratings on linear and circular stimulus dimensions. QUEST+ provides a general method to accelerate data collection in many areas of cognitive and perceptual science.**

The experimenter assumes a parametric form of the psychometric function … After each response, values of the parameters are computed that maximize the probability of the set of responses that have been obtained, given the set of stimuli that have been presented. These values are then used to determine the placement of the next stimulus.

*failure*and

*success*,

*incorrect*and

*correct*,

*no*and

*yes*, or another congruent pair. But the Bayesian framework is not limited to that case. Below we show how other psychometric procedures, such as categorization, with more than two possible outcomes can be handled by QUEST+.

**x**, the vector of psychometric function parameters is written

**s**, and the outcome is written

*r*. The psychometric function itself is written

**Ψ(x,s)**. In contrast to previous treatments, here the psychometric function always returns a vector of probabilities, corresponding to the several possible outcomes

*r*. For a given

**x**and

**s**, those probabilities must sum to one. Alternatively, we can express the psychometric function as the conditional probability of outcome

*r*given

**x**and

**s**or

*p*(

*r|*

**x**

*,*

**s**). The stimulus value on trial

*k*is a vector

**x**

*, and the complete sequence of stimulus values is*

_{k}**X**

*= {*

_{k}**x**

_{1},

**x**

_{2},…,

**x**

*}. The corresponding sequence of responses is*

_{k}**r**

*= {*

_{k}*r*,

_{1}*r*, …,

_{2}*r*}. Although it is arbitrary, we generally consider the several possible outcomes to be identified by successive integers starting at one (e.g.,

_{k}*r*= 1 -> failure,

*r*= 2 -> success). We adopt the convention of using uppercase

*P*to indicate probability density functions and lowercase

*p*to indicate likelihoods or conditional probabilities. We also use the prime symbol to indicate un-normalized density functions; the corresponding symbol without a prime is understood to have been normalized by dividing by its integral.

*P*(

**s**). From Bayes' theorem, we can write the posterior probability of the parameters after trial

*k*as the product of the prior density and the likelihood of the data Thus, the posterior density after trial k + 1 is given by In words, we multiply the posterior density on trial

*k*by the likelihood of trial

*k +*1.

**c**and

**threshold**are contrasts in dB, and

**slope**,

**guess**, and

**lapse**are the traditional parameters of the Weibull function (Treutwein, 1995; Watson & Pelli, 1983).

**Outcome[x_List]**, where

**x**is a list of stimulus parameters and the function must return one of the possible outcomes. The outcomes are identified by positive integers starting at one. The number of possible outcomes

*n*is determined by the psychometric function. For the Weibull function above,

*n =*2.

**x**, the name of a psychometric function

**psi**, and a list of psychometric parameter definitions

**s**. This function makes use of the multinomial distribution, which is a generalization of the binomial distribution. The psychometric function returns the

*n*probabilities of the

*n*possible responses, and the SimulatedObserver function randomly returns a number between one and

*n*with the specified probabilities.

**x**. In our implementation, we do this by specifying the actual samples in each dimension. We use the Mathematica option construct to specify these quantities. Here is an example in which we indicate that the single stimulus parameter should take values from −40 to zero in steps of one. Note that the structure is a list of lists because QUEST+ allows more than one stimulus parameter.

**s**. In this case, we vary two parameters: threshold and slope. But the Weibull psychometric function defined above has four parameters. Here we handle that by specifying single values (lists of length one) for the fixed parameters. So the specification is a list of four lists. The first is for threshold and the second for slope, and the remaining two are fixed values for guess and lapse.

**QpPriors**option. The default setting is

options = Sequence[

QpOutcome -> (SimulatedObserver[#, QpPFWeibull, {-20, 3, .5, .02}] & )

, QpPNames -> {"Threshold", "Slope", "Guess", "Lapse"}

, QpXNames -> {"Contrast (dB)"}

, QpXSamples -> {Range[-40, 0]}

, QpPSamples -> {Range[-40, 0], Range[2, 5], {.5}, {0.02}}

, QpPsi -> QpPFWeibull

, QpPrior -> "Constant"];

**QpPNames**and

**QpXNames**, as they are used only for subsequent analysis and plotting or are defaults, such as

**QpPrior**, but we show them here for completeness.

result = QuestPlus[32, options]{{-20,3,0.5,0.02},{{{-18},2},{{-22},1},{{-12},2},{{-13},2},{{-15},2},{{-16},2},{{- 17},2},{{-18},2},{{-19},2},{{- 21},2},{{-24},2},{{-25},2},{{-26},2},{{-27},1},{{- 21},2},{{-22},2},{{-23},1},{{-20},2},{{-20},2},{{-20},2},{{-21},1},{{-19},2},{{- 19},2},{{-20},2},{{-20},2},{{-20},1},{{-19},2},{{-19},2},{{-19},2},{{-19},2},{{- 19},2},{{-20},2},{{-20},2},{{-20},2},{{-20},2},{{-20},2},{{-21},2},{{-21},1},{{- 20},2},{{-20},2},{{-20},1},{{-19},1},{{-18},2},{{-18},2},{{-19},2},{{-19},2},{{- 19},1},{{-18},2},{{-18},1},{{-17},2},{{-18},2},{{-18},2},{{-18},1},{{-17},2},{{- 17},2},{{-17},2},{{-17},2},{{-17},2},{{-17},2},{{-17},2},{{-18},2},{{-18},2},{{- 18},2},{{-18},2}}}

**QuestPlus**function above assumes that the user provides only a single additional function: the function that presents the stimulus to the observer. In that sense,

**QuestPlus**is the “master” program. Users may prefer to write their own master program and invoke QUEST+ functions as subordinate routines. This is useful when a number of conditions are to be interleaved, when stimulus creation and presentation are complex, or when the user would like to deviate from the stimulus parameters for the next trial recommended by QUEST+. We provide this capability in the form of four functions:

**QpInitialize**,

**QpQuery**,

**QpUpdate**, and

**QpRun**. The last is only an example of how to invoke the first three, and we show it here in its entirety.

QpRun[ntrials_, options___Rule] :=Block[{outcomefunction, outcome, structure, psamples},

{outcomefunction, psamples} = {QpOutcome, QpPSamples} /. {options} /. Options[QuestPlus];

structure = QpInitialize[options];

Do[(

structure=QpQuery[structure];

outcome=outcomefunction[structure[[5]]];

structure=QpUpdate[structure, outcome];)

,{ntrials}];

{MapThread[#1[[#2]]&,{psamples, ListArgMax[structure[[2]]]}],structure[[1]]}]

**QuestPlus**,

**QpRun**returns the quantized parameter estimates and the data.

**QpInitialize**creates a data structure that is subsequently queried and updated on each trial. When several conditions are interleaved,

**QpInitialize**would be called once for each to create an array of structures.

**QpQuery**determines the recommended next stimulus vector, and

**QpUpdate**uses the previous outcome to update the data structure.

**QpFit**fits a psychometric function to the result. This will typically be the psychometric function used in the QUEST+ procedure and specified in the options. Although QUEST+ effectively finds the ML fit on every trial,

**QpFit**is not constrained to the quantized parameter space. But because it takes advantage of the same estimation logic, it operates on any number of stimulus dimensions, psychometric parameters, and outcomes.

**QpPlot**,

**QpPlot2D**,

**QpPlot3D**, and

**QpPlotCategories**to show the results of one, two, and three stimulus dimensions and of categorization (rating) experiments. These functions first fit the data using

**QpFit**and then plot data and best-fitting psychometric function. Examples of these plots are shown below.

Options[QuestPlus] =

{ QpOutcome -> (SimulatedObserver[#, QpPFWeibull, {-20, 3.5, .5, .02}] &)

, QpPNames -> {"Threshold (dB)", "Slope", "Guess", "Lapse"}

, QpXNames -> {"Contrast (dB)"}

, QpXSamples -> {Range[-40, 0]}

, QpPSamples -> {Range[-40, 0], {3.5}, {.5}, {0.02}}

, QpPsi -> QpPFWeibull

, QpPrior -> "Constant"

, Verbose -> False};

result = QuestPlus[32]{{-20, 3.5, 0.5, 0.02}, {{{-18}, 2}, {{-22}, 2}, {{-25}, 2}, {{-28}, 2}, {{-30}, 1}, {{-22}, 1}, {{-13}, 2}, {{-15}, 2}, {{-16}, 2}, {{-18}, 2}, {{-19}, 2}, {{-20}, 2}, {{-21}, 2}, {{-22}, 2}, {{-23}, 1}, {{-19}, 2}, {{-20}, 2}, {{-20}, 1}, {{-18}, 2}, {{-18}, 2}, {{- 19}, 1}, {{-17}, 2}, {{-17}, 2}, {{-18}, 2}, {{-18}, 2}, {{-18}, 2}, {{-19}, 2}, {{-19}, 2}, {{-19}, 2}, {{-19}, 2}, {{-19}, 2}, {{-19}, 2}}}

options = Sequence[

QpOutcome -> (SimulatedObserver[#, QpPFWeibull, {-20, 3, .5, .02}] & )

, QpXNames -> {"Contrast (dB)"}

, QpPNames -> {"Threshold (dB)", "Slope", "Guess", "Lapse"}

, QpXSamples -> {Range[-40, 0]}

, QpPSamples -> {Range[-40, 0], Range[2, 5], {.5}, Range[0, .04, .01]}

, QpPsi -> QpPFWeibull];

*c*is the contrast threshold in dB,

*f*is spatial frequency in cycles per degree,

*w*is temporal frequency in Hz,

*t*is the minimum threshold, and

*c*,

_{0}*c*, and

_{f}*c*are the linear coefficients. If we set

_{w}*w*= 0, we obtain the spatial contrast threshold function (the inverse of the contrast sensitivity function) as shown in Figure 5A.

*c*,

_{0}*c*, and

_{f}*t*; and finally, the psychometric function whose parameters are to be estimated.

options = Sequence[

QpOutcome -> (SimulatedObserver[#, QpPFSCSF, {-35, -50, 1.2}] & )

, QpPNames -> {"t", "c_{0}", "c_{f}"}

, QpXNames -> {"Spatial frequency (c/d)", "Contrast (dB)"}

, QpXSamples -> {Range[0, 40, 2], Range[-50, 0, 2]}

, QpPSamples -> {Range[-30, -50, -2], Range[-40, -60, -2], Range[.8, 1.6, .2]}

, QpPsi -> QpPFSCSF];

result = QuestPlus[32, options]{{-32,-56,1.4},{{{40,0},2},{{40,-4},1},{{34,0},2},{{36,0},2},{{0,-38},2},{{38,0},2},{{0,- 40},1},{{38,0},2},{{18,- 26},2},{{18,-26},2},{{40,0},2},{{0,-36},2},{{0,- 36},1},{{40,0},2},{{20,-26},2},{{20,-26},2},{{40,-2},2},{{22,-26},2},{{40,- 6},1},{{20,-26},1},{{40,0},2},{{18,-26},2},{{40,0},2},{{18,-26},2},{{0,-32},2},{{0,- 32},2},{{0,-34},2},{{0,-34},1},{{40,0},1},{{18,-26},2},{{38,0},2},{{18,-26},2}}}

*c*.

_{f}

options=Sequence[

QpOutcome-> (SimulatedObserver[#, QpPFSTCSF,{-40,-50,1.2,1.0}]& )

, QpPNames->{"t”,"c_{0}”,"c_{f}”,"c_{w}"}

, QpXNames->{"Contrast (dB)”,"Spatial frequency (c/d)”,"Temporal frequency (Hz)"}

, QpXSamples->{Range[-40,0,5], Range[0,40,5], Range[0,40,5]}

, QpPSamples->{Range[-30,-50,-5], Range[-40,-60,-5], Range[.8,1.6,.2], Range[.8,1.6,.2]}

, QpPsi-> QpPFSTCSF];

**QpPlot3D**. In order to see the placement of the trials, we show the three-dimensional stimulus space as a series of plots, one for each of the contrast samples. Each location tested is plotted as a disc with a size reflecting the number of trials at that location, and the angular portion in green represents the percentage correct. We superimpose these points on gray-level representations of the psychometric function of the simulated observer. This representation of the results is shown in Figure 9. QUEST+ again concentrates the trials at the critical locations in the three-dimensional stimulus space. In 64 trials, of the 729 possible stimulus locations, only 20 are visited. In a subsequent simulation of 256 trials, only 24 were visited.

*x*and perceptual values is given by a scaling function Θ(

*x*).

*SD*of 1. The internal scale is effectively measured in units of the standard deviation, so we can also consider it to express units of just noticeable difference or

*d′*. Consider an experiment consisting of pair comparisons: asking the observer which is greater of two stimuli differing along the physical dimension. According to the Thurstone model, the probability of a correct judgment in a pair comparison is a function only of the distance between the sensory magnitudes induced by the two intensities of the pair where Φ is the cumulative Normal distribution.

*x*

_{1}and

*x*

_{2}. However, this is typically a daunting task because, if there are

*n*discrete values along the

*x*dimension and we need

*k*trials at each to get reliable results, we will need

*k n*

^{2}trials. And many trials will be uninformative because they involve separations that are much too easy to discriminate or much too hard. Previously, we designed an adaptive method to deal efficiently with this problem (Watson & Kreslake, 2001). Here we show that it fits within the QUEST+ methodology.

*m*, the maximum scale value,

*t*, a threshold value of

*x*below which the response is zero, and

*q*, a power. This model is widely used in psychophysical scaling theory (Engen, 1972; Falmagne, 1986). An example is illustrated in Figure 10. The function is zero until the threshold is reached, and then it rises as a power function to the maximum at

*x =*1.

*x*

_{2}= 1.

*x*

_{2}is larger than

*x*

_{1}, but because both are values from the same dimension, we get the same information from a trial at {

*x*

_{1},

*x*

_{2}} as we get from {

*x*

_{2},

*x*

_{1}}. Because we are only interested in deriving the parameters, it doesn't matter to us whether Quest+ chooses one or the other. Likewise, we could process the data in such a way that on all trials in which

*x*

_{1}is greater than

*x*

_{2}, we exchange the stimulus values and exchange the responses.

*δ*(arcmin) is given by where

*σ*is the reference blur (arcmin),

*β*(arcmin) is an intrinsic “neural” blur, and

*ω*is a parameter greater than one. An example of this function is shown in Figure 12 on log–log coordinates.

options = Sequence[

QpOutcome -> (SimulatedObserver[#, QpPFBlur, {1.2, 1}] & )

, QpPNames -> {"omega”, “beta"}

, QpXNames -> {"Reference blur (log arcmin)”, “Threshold blur (log arcmin)"}

, QpXSamples -> {Range[-1, 1.5, .25], Range[-1, 1, .25]}

, QpPSamples -> {Range[1.05, 1.5, .05], Range[.5, 2, .1]}

, QpPsi -> QpPFBlur];

(result = QuestPlus[32,options])//Short{{1.2,1.},{{{-0.25,-0.25},2},<<31>>}}

*SD*of 0.5 and category widths of {1, 2, 3}, meaning that the criteria are at one, three, and six. The psychometric function returns four probabilities, corresponding to the four responses. At any stimulus value, of course, the four curves must sum to one.

options = Sequence[

QpOutcome -> (SimulatedObserver[#, QpPFRating, {.5, 1, 2, 3}] & )

, QpXNames -> {"Magnitude"}

, QpPNames -> {"SD", "Width1", "Width2", "Width3"}

, QpXSamples -> {Range[0, 9, .25]}

, QpPSamples -> {Range[.05, .5, .05], Range[.5, 4, .5], Range[.5, 4, .5], Range[.5, 4, .5]}

, QpPsi -> QpPFRating];

**QpPlotCategories**. To allow us to show all four responses, each is shown with a different quarter circle symbol so that they do not overlap. As in Figures 1 through 4, the color saturation reflects the number of trials. The trial placement reflects the strategy of entropy minimization: Few trials are placed in the center of any category; rather they are placed near the boundaries (the criteria).

*n*categories on a circular dimension extending from 0 to 2

*π*. We assume the internal response follows a von Mises distribution with parameters

*mean*and

*concentration*(Fisher, 1995). Concentration is similar to inverse variance, and when it is zero, the density is constant at 1/(2

*π*). We assume the mean is given by the stimulus value. We further assume that the observer maintains

*n*− 1 criteria that separate the circular dimension into

*n*regions associated with the

*n*possible responses.

*n*probabilities. The list of parameters is {concentration, criterion

_{1}, width

_{1}, …, width

_{n}_{−}

_{1}}. Only

*n*− 1 widths are specified because the

*n*widths must total 2

*π*. We can then configure QUEST+ with the following options:

options = Sequence[

QpOutcome -> (SimulatedObserver[#, QpPFCircular, {4, Pi/3, 2 Pi/3, 2 Pi/3}] & )

, QpNames -> {"Concentration", "c1", "w1", "w2"}

, QpXNames -> {"Direction"}

, QpXsamples -> {Range[0, 17] Pi/9}

, QpPsamples -> {Range[1, 8], Range[0, 17] Pi/9, {2 Pi/3}, {2 Pi/3}}

, QpPsi-> QpPFCircular];

*π*/3, and equal widths of 2

*π*/3. Because two widths are specified, there will be three possible outcomes. In the list of parameter samples, note that for each of the two width parameters we give only a single value. This means the parameter will not be estimated but fixed at that value. This is a convenient way of specifying a fixed parameter.

*Define the default options.*

Options[QuestPlus] = {

QpOutcome -> (SimulatedObserver[#, QpPFWeibull, {-20, 3.5, .5, .02}] &),

QpPNames -> {"Threshold (dB)”, "Slope", "Guess", "Lapse"},

QpXNames -> {"Contrast (dB)"},

QpXSamples -> {Range[-40, 0]},

QpPSamples -> {Range[-40, 0], {3.5}, {.5}, {0.02}},

QpPsi -> QpPFWeibull,

QpPrior -> "Constant"};

*The main function.*

*Define the local variables.*

*Set the values of the variables that can be set by options.*

*Initialize the list of trials to an empty list.*

*Compute the number of dimensions of the stimulus array.*

*Compute the number of dimensions of the parameter array.*

*Compute the number of dimensions of the likelihood array.*

*Compute the depth in the likelihood array at which parameter array begins.*

*If necessary, define a constant prior. It is a uniform array with the same dimensions as the parameter array.*

*Compute the likelihood array (conditional probabilities of each outcome at each stimulus value for each parameter value).*

*Set the posterior array equal to the prior array.*

*For each trial…*

*Compute the product of likelihood and current posterior array at each outcome and stimulus location.*

*Compute the total probability of each outcome at each stimulus location.*

*Compute posterior PDF for each stimulus location and outcome.*

*Compute the expected entropy for each stimulus location.*

*Find the index of the stimulus with the smallest expected entropy.*

*Set the next stimulus location to the location of the smallest expected entropy.*

*Do a trial of a simulated or real observer.*

*Update the list of trials.*

*Set the new posterior to one of the previously computed posteriors, based on the outcome and stimulus location.*

*Repeat until the number of trials is complete.*

*Estimate the optimal quantized Bayesian parameters by the location of the maximum in the final posterior.*

*M*, where

*N*is an operator giving the number of samples in an array. This is just the product of the lengths of all the lists of responses, of stimulus samples, and parameter samples. Clearly, this number can become large rather quickly if care is not taken. However, on fast modern computers, this direct approach is quite practical for many problems. We have performed timing tests on the 10 examples discussed above, and the results show that the time per trial is constant for each experiment. The timings below were conducted in the relatively slow Mathematica language (Version 10.4) on a modest laptop (MacBook Pro, 2.5 GHz Intel Core i7 processor, 16 GB 1600 MHz DDR3 memory). We conducted simulations with 32, 64, and 128 trials and fit the resulting timings with a linear function of trials. In each case, the fit was nearly perfect. The slope of each fit is a measure of seconds required to compute the stimulus location for the next trial. These values are shown as s/trial in Table 1. These values are always less than 1 s/trial for the examples tested. We also divide s/trial by the number of samples in the likelihood array (

*M*), and that is shown in Table 1 as

*μ*s/trial/sample.

*, 10996.*

*Nature Communications 7**, 15 (9): 5, 1–20, doi:10.1167/15.9.5. [PubMed] [Article]*

*Journal of Vision**, 39, 151–153.*

*Perception & Psychophysics**(pp. 47–86). New York: Holt, Rinehart and Winston.*

*Woodworth and Schlosberg's Experimental Psychology**(pp. 1.3–1.66). New York: Wiley.*

*Handbook of Perception and Human Performance**. Cambridge, UK: Cambridge University Press.*

*Statistical analysis of circular data**, 11 (6), 1710–1719.*

*Journal of the Optical Society A**, 44, 370.*

*Journal of the Acoustical Society of America**, 34 (7), 885–912.*

*Vision Research**, 37 (12), 1595–1604.*

*Vision Research**, 39 (16), 2729–2737.*

*Vision Research**, 50 (4), 369–389.*

*Journal of Mathematical Psychology**, 70 (12), 1458–1471.*

*Journal of the Optical Society of America**, 9 (8): 696, doi:10.1167/9.8.696. [Abstract]*

*Journal of Vision**, 46 (19), 3160–3176.*

*Vision Research**, 10 (3): 17, 1–21, doi:10.1167/10.3.17. [PubMed] [Article]*

*Journal of Vision**, 60 (3), 382–389.*

*Journal of the Optical Society of America**, 28( 3, Suppl.), 366.*

*Investigative Ophthalmology & Visual Science**, 16 (3), 647–653.*

*Journal of the Optical Society of America A**, 57 (3), 773–786.*

*Journal of Neurophysiology**, 41 (4), 431–439.*

*Journal of Mathematical Psychology**, 13 (10): 9, 1–8, doi:10.1167/13.10.9. [PubMed] [Article]*

*Journal of Vision**. Chicago: University of Chicago Press.*

*The measurement of values**, 35 (17), 2503–2522.*

*Vision Research**, 23 (5), 483–515.*

*Seeing and Perceiving**, 11 (5): 10, 1–23, doi:10.1167/11.5.10. [PubMed] [Article]*

*Journal of Vision**, 2016 (16), 1–6, doi:10.2352/ISSN.2470-1173.2016.16.HVEI-102.*

*Electronic Imaging**, 4299, 79–89.*

*Proceedings of the SPIE**, 14, 6–7.*

*Applied Vision Assoc. Newsletter**, 33 (2), 113–120.*

*Perception & Psychophysics**, 453 (7192): 233–235.*

*Nature*