^{1}

^{1}

^{2}

^{2}

^{3}

^{*}

We consider the response of a test subject upon a skin area being heated with an electromagnetic wave or a contact surface. When the specifications of the electromagnetic beam are fixed, the stimulus is solely described by the heating duration. The binary response of a subject, escape or no escape, is determined by the stimulus and a subjective threshold that varies among test realizations. We study four methods for inferring the median subjective threshold in psychophysical experiments: 1) sample median, 2) maximum likelihood estimation (MLE) with 2 variables, 3) MLE with 1 variable, and 4) adaptive Bayesian method. While methods 1 - 3 require samples of time to escape measured in the method of limits, method 4 utilizes binary outcomes observed in the method of constant stimuli. We find that a) the adaptive Bayesian method converges and is as efficient as the sample median even when the assumed model distribution is incorrect; b) this robust convergence is lost if we infer the mean instead of the median; c) for the optimal performance in an uncertain situation, it is best to use a wide model distribution; d) the predicted error from the posterior standard deviation is unreliable, dominated by the assumed model distribution.

In psychophysical experiments, it is always desirable to understand the mechanisms that give rise to the observed behavior [

We consider the response of a test subject upon a skin area being heated with a contact surface [

Due to the biovariability and other uncertainties, the escape time of the test subject will be different from one trial to another. We model the subject response as solely determined by the stimulus (the exposure duration) and a subjective threshold that varies among test subjects and varies among different trial realizations for a given subject. The outside stimulus and the internal threshold are connected via the pain signal generated at the heat-sensitive nociceptors activated by the skin temperature increase [

Let Z be the subjective threshold. We set t = 0 as the start of heating exposure. By time t = u , the cumulative stimulus (the exposure duration) is u and the subject response is determined by stimulus u and subjective threshold Z.

Subject response ( u ) = { no escape if u < Z escape if u > Z (1)

In a particular trial realization, Z is fixed but unknown. In (1), we write the subject response as a function of stimulus u to emphasize that in a trial, as the exposure duration u increases, the escape response will eventually occur. Over many subjects and over many trial realizations, mathematically Z is a random variable, reflecting the biovariability and other uncertainties. By definition, Z is the minimum exposure duration needed to induce the escape response for a particular subject in a particular trial. Since Z is naturally constrained by Z > 0 , we use a Weibull distribution to model Z. The actual distribution of Z may be different. We will investigate the consequence of this discrepancy on the inference accuracy. The probability density function (PDF) of the Weibull distribution is

ρ Z ( z ) = exp ( − z β α β ) β z β − 1 α β , z > 0 (2)

which is described by two parameters: shape parameter β and scale parameter α . The mean, median and standard deviation of U have the expressions

E ( Z ) = α Γ ( 1 + 1 β )

z ( m ) ≡ med ( Z ) = α ( ln 2 ) 1 β (3)

std ( Z ) = α Γ ( 1 + 2 β ) − Γ ( 1 + 1 β ) 2

Here med ( Z ) denotes the median of Z, and Γ ( ) is the gamma function.

Two types of psychophysical experiments may be conducted using, respectively, the method of limits (ML) [

Let z be the sample value of the subjective threshold in a particular trial. Under the condition that the test subject is unaware of the test type used in the trial, it is reasonable to expect that threshold z is independent of the test type. In a test using the method of constant stimuli, let u be the prescribed stimulus (exposure duration). For a random trial with stimulus u, under the assumption that Z has a Weibull distribution, the probability of escape is given by the Weibull psychometric function denoted by P ( u ) , which is the cumulative distribution function (CDF) of the Weibull distribution.

Pr ( escape ) = Pr ( Z < u ) = ∫ 0 u ρ Z ( z ) d z = 1 − exp ( − u β α β ) ≡ P ( u ) (4)

The median subjective threshold z ( m ) ≡ med ( Z ) satisfies P ( z ( m ) ) = 50 % . In the context of psychometric function P ( u ) , u = z ( m ) is the point of subjective equality (PSE), the stimulus level at which it is equally likely that the heat induced is above or below the subjective pain tolerance. If we view P ( u ) as a dose-response relation, u = z ( m ) is the half maximum response dose. We write the psychometric function P ( u ) in terms of z ( m ) and β .

P ( u ) = 1 − exp ( − ( ln 2 ) u β ( z ( m ) ) β ) (5)

The goal of this study is to estimate z ( m ) from data/observations. We choose to infer z ( m ) = med ( Z ) instead of E ( Z ) for several reasons:

· z ( m ) is the point of subjective equality (PSE) in the psychometric function.

· If we measure samples of Z using the method of limits, the sample median is less susceptible to outliers, which may occur in real experiments.

· The median is preserved in any monotonic transformation, which allows us to work with random variable Y = log 10 Z to infer med ( Y ) in the adaptive Bayesian method. Notice that while Z is constrained by Z > 0 , variable Y is unconstrained.

· In the adaptive Bayesian method, the convergence of the inferred median is independent of whether the assumed model is correct. In our simulation analysis, we will demonstrate that even when the assumed model deviates from the actual psychometric function, the inferred median still converges to the correct value. If we choose E ( Z ) as the inference variable, however, this convergence property is lost.

We discuss four inference methods. Each method is based on different assumptions and thus is applicable in different situations.

This method does not assume any particular distribution form of Z. It is applicable for estimating the median of any random variable. We measure independent samples of subjective threshold Z in psychophysical tests using the method of limits. The data set is

Z ( d a t a ) = { z j , j = 1 , 2 , ⋯ , N } (6)

We use z ^ ( m ,1 ) to denote the inference result of method 1 for med ( Z ) .

z ^ ( m ,1 ) = samplemedianof { z j , j = 1,2, ⋯ , N } (7)

Like method 1 above, the maximum likelihood estimation (MLE) requires data set Z ( data ) , consisting of measured independent samples of Z. MLE (2 variables) assumes that the subjective threshold Z has the Weibull distribution given in (2), with shape parameter β and scale parameter α both unknown. MLE (2 variables) is applicable only when Z has the Weibull distribution. MLE (2 variables), we infer ( α , β ) simultaneously by maximizing the likelihood and then calculate z ( m ) from ( α , β ) . For mathematical convenience, let η ≡ α β . In terms of ( η , β ) , the likelihood function and the log-likelihood have the expressions:

L ( η , β ; Z ( data ) ) = ∏ j = 1 N ρ Z ( z j ) = ∏ j = 1 N exp ( − z j β η ) β z j β − 1 η

l ( η , β ; Z ( data ) ) = ln L ( η , β ; Z ( data ) ) = − 1 η ∑ j = 1 N z j β − N ln η + ∑ j = 1 N ln ( β z j β − 1 ) (8)

We maximize the log-likelihood l ( η , β ; Z ( data ) ) to estimate ( η ^ , β ^ ) .

( η ^ , β ^ ) = arg max ( η , β ) l ( η , β ; Z ( data ) ) (9)

Once η ^ and β ^ are determined, the inference result of method 2 for med ( Z ) , denoted by z ^ ( m ,2 ) , is calculated using Equation (3) and α ^ = η ^ 1 β ^ .

z ^ ( m ,2 ) = ( η ^ ln 2 ) 1 β ^ (10)

Like methods 1 and 2 above, method 3 requires data set Z ( data ) , consisting of measured independent samples of Z. This method assumes 1) the subjective threshold Z has the Weibull distribution of the form (2) and 2) the shape parameter β is known, which leaves the scale parameter α as the only unknown. MLE (1 variable) is applicable only when Z has the Weibull distribution and the true value of β is given. We write the log-likelihood in terms of η ≡ α β .

l ( η ; Z ( data ) ) = ln L ( η ; Z ( data ) ) = − 1 η ∑ j = 1 N z j β − N ln η + ∑ j = 1 N ln ( β z j β − 1 ) (11)

The log-likelihood has the derivative

d d η l ( η ; Z ( data ) ) = 1 η 2 ∑ j = 1 N z j β − N η .

The maximization problem of l ( η ; Z ( data ) ) has an analytical solution.

η ^ = arg max η l ( η ; Z ( data ) ) = 1 N ∑ j = 1 N z j β (12)

The inference result of method 3 for med ( Z ) , denoted by z ^ ( m ,3 ) , is calculated using (3)

z ^ ( m , 3 ) = η ^ 1 β = ( 1 N ∑ j = 1 N z j β ) 1 β (13)

Unlike methods 1-3 above, the adaptive Bayesian method does not require measured independent samples of Z. Rather it utilizes binary outcomes observed in tests using the method of constant stimuli (MCS). Suppose a sequence of n trials of the MCS type have already been completed with prescribed stimuli (exposure durations) { u j , j = 1 , 2 , ⋯ , n } . The observed binary outcome I j of each trial indicates either z j > u j (if no escape occurs) or z j < u j (if escape occurs).

I j = { 0 ( no escape ) if z j > u j 1 ( escape ) if z j < u j (14)

Here z j is the (unknown) value of subjective threshold in trial j. In trials of the MCS type, values of { z j } are not fully revealed. Instead, only z j > u j or z j < u j is observed.

Following the approach used in the Quest package [

D ( n ) = { ( x j , I j ) , j = 1 , 2 , ⋯ , n } (15)

To infer z ( m ) ≡ med ( Z ) in a Bayesian framework, we treat y ( m ) ≡ log 10 z ( m ) as a random variable and use a normal distribution as the prior.

Priorof y ( m ) ∼ N ( μ y , σ y 2 ) (16)

For a trial with stimulus x (exposure duration u = 10 x ), under the assumption that the subjective threshold Z has a Weibull distribution, the probability of escape is described by the psychometric function P ( ⋅ ) in (5)

Pr ( I = 1 ) = Pr ( Z < 10 x ) = P ( 10 x ) = 1 − exp ( − 10 β ( x − y ( m ) ) + δ ) ≡ F ( x ) , δ ≡ log 10 ( ln 2 )

The model psychometric function F ( x ) has parameters β and y ( m ) = log 10 z ( m ) . In experiments, however, the probability of escape is governed by the actual psychometric function F actual ( x ) , which may be different from F ( x ) .

In the adaptive Bayesian method, we assume β is given and our goal is to infer y ( m ) . For that purpose, we need to view the probability of escape as a function of y ( m ) . We write F ( x ) in terms of ( x − y ( m ) ) and β as.

F ( x ) = ϕ ( x − y ( m ) ; β ) (17)

where ϕ ( s ; β ) is parameter-free function defined as

ϕ ( s ; β ) ≡ 1 − exp ( − 10 β s + δ ) , δ ≡ log 10 ( ln 2 ) (18)

To facilitate the inference of y ( m ) , we shift y ( m ) by its prior mean and consider ξ ≡ y ( m ) − μ y . The prior distribution of ξ is

ρ ξ ( pri ) ( ξ ) ∼ exp ( − ξ 2 2 σ y 2 ) (19)

(19) is our prior knowledge about ξ . When data set D ( n ) from n trials becomes available, the posterior distribution of ξ is calculated according to the Bayes theorem.

ρ ξ ( post ) ( ξ | D ( n ) ) ∼ ρ ξ ( pri ) ( ξ ) ∏ j = 1 n B ( ϕ ( x j − μ y − ξ ; β ) , I j ) (20)

where B ( p , k ) ≡ p k ( 1 − p ) 1 − k is the Bernoulli distribution. A key component of the adaptive Bayesian method is to select the stimulus for the next trial based on the current posterior of ξ . Supposen trials have been conducted and the data set D ( n ) is used to calculate the posterior ρ ξ ( post ) ( ξ | D ( n ) ) in (20). Let x n + 1 denote the stimulus for the next trial. The optimal choice of x n + 1 is the point of subjective equality (PSE) of the actual psychometric function where F actual ( x n + 1 ) = 50 % . The PSE of the actual psychometric function is unknown. It is the inference variable. We use the best available estimate of PSE as the next stimulus. We set stimulus x n + 1 to the posterior median of y ( m ) = μ y + ξ , which reflects the prior of y ( m ) and the observed outcomes of the preceding n trials.

x n + 1 = μ y + med ( ρ ξ ( post ) ( ξ | D ( n ) ) ) (21)

Once stimulus x n + 1 is prescribed, trial ( n + 1 ) is carried out and the observed binary outcome I n + 1 is used to update the posterior distribution of ξ .

ρ ξ ( post ) ( ξ | D ( n + 1 ) ) ∼ ρ ξ ( post ) ( ξ | D ( n ) ) B ( ϕ ( x n + 1 − μ y − ξ ; β ) , I n + 1 ) (22)

It is important to emphasize that the probability of escape in real experiments is governed by Pr ( I n + 1 = 1 ) = F actual ( x n + 1 ) , independent of the assumed model distribution.

We repeat the process of adaptively selecting stimulus x based on the current posterior distribution, carrying out the next trial and using the observed outcome to update the posterior distribution. The iterative process of (21)-(22) is repeated either for a specified number of trials or until the standard deviation of posterior ξ is less than a specified tolerance. When the adaptive Bayesian iteration is completed, the inference result of method 4 for med ( Z ) , denoted by z ^ ( m ,4 ) , is calculated as the posterior median of z ( m ) = 10 y ( m )

z ^ ( m ,4 ) = 10 y ^ ( m ,4 ) , y ^ ( m ,4 ) = μ y + med ( ρ ξ ( post ) ( ξ | D ( N ) ) ) (23)

We use Monte Carlo simulations to evaluate the inference accuracy for each of the four methods described in the previous section.

We define the relative error as the root-mean-square (RMS) relative difference between the inference result and the true value. For each inference method, we carry out M rounds of Monte Carlo simulations. For methods 1 - 3, each round of simulations involves generating a data set of N samples and calculating an inference result from each simulated data set. For method 4, each round of simulations consists of repeating for N trials the Bayesian process of adaptively selecting the stimulus based on the posterior distribution, carrying out the trial and updating the posterior distribution according to the observed outcome. Let

· z ( m , e ) be the true value of med ( Z ) and

· z ^ ( m , k , i ) be the inferred med ( Z ) , using method k in round i of simulations.

From the Monte Carlo simulations, we calculate the relative error of method k as

err ( k ) = 1 M ∑ i = 1 M ( z ^ ( m , k , i ) − z ( m , e ) z ( m , e ) ) 2 (24)

We use M = 10 6 rounds of Monte Carlo simulations in the calculation of each relative error (each entry in

Before we select values of ( α , β ) , we point out that random variable Z, the exact and inferred (mean, median, std) of Z, are all proportional to scale parameter α . In particular, both z ^ ( m , k , i ) and z ( m , e ) are proportional to α . Thus, the relative error is independent of α .

In simulations of this section, we use β = 2.5 and 5. For the convenience of setting up simulations, we set α = ( ln 2 ) − 1 β to make med ( Z ) = 1 .

In the adaptive Bayesian method, an estimate of z ( m ) is obtained via inferring y ( m ) ≡ log 10 z ( m ) . The prior for y ( m ) is a normal distribution: y ( m ) ∼ N ( μ y , σ y 2 ) . The corresponding prior of z ( m ) is a log-normal distribution.

ρ ( z ( m ) = s ) = 1 q 0 s 2 π σ y 2 exp ( − ( log 10 s − μ y ) 2 2 σ y 2 ) , q 0 = ln 10 (25)

The prior (mean, median and standard deviation) of z ( m ) are related to

( μ y , σ y ) by

E ( z ( m ) ) = 10 μ y exp ( q 0 2 σ y 2 ) , q 0 = ln 10

μ z ≡ med ( z ( m ) ) = 10 μ y

σ z ≡ std ( z ( m ) ) = 10 μ y exp ( q 0 2 σ y 2 ) [ exp ( q 0 2 σ y 2 ) − 1 ] (26)

Suppose we know the prior median and standard deviation of z ( m ) from existing data of measurements on the subjective threshold Z. Given ( μ z , σ z ) , we can build the normal prior of y ( m ) by solving for ( μ y , σ y ) in (26)

μ y = log 10 μ z

σ y = 1 q 0 ln ( 1 2 + 1 2 1 + 4 ( σ z μ z ) 2 ) , q 0 = ln 10 (27)

For small ( σ z / μ z ) , the Taylor expansion of the RHS yields

σ y = 1 q 0 ( σ z μ z ) [ 1 − 3 4 ( σ z μ z ) 2 + ⋯ ] ≈ 1 ln 10 ( σ z μ z ) (28)

The prior relative uncertainty of z ( m ) is described by the ratio σ z / μ z . From the given σ z / μ z , the prior standard deviation of y ( m ) is calculated using Equation (27). In Monte Carlo simulations, the prior mean of y ( m ) is set to a random variable

μ y ∼ N ( y ( m , e ) , σ y 2 )

where y ( m , e ) = log 10 z ( m , e ) is the exact median of y ( m ) = log 10 z ( m ) . This randomness of μ y in simulations is both realistic and necessary. In real applications, the exact value of y ( m , e ) is unknown. In the adaptive Bayesian method, the inference result for y ( m ) is given by the median of the posterior y ( m ) . If we set μ y = y ( m , e ) , then before any Bayesian updating, the prior distribution already gives the exact solution, which is unreasonable. Setting μ y ∼ N ( y ( m , e ) , σ y 2 ) is consistent with that the prior error of y ( m ) is σ y .

In all simulations, data sets are generated using the actual distribution. In Section 4, the actual distribution and the assumed model distribution have the same Weibull form and the same shape parameter β . The median of the actual distribution, z ( m , e ) , is set to 1 while the median of the model distribution, z ( m ) , is the inference variable. Not all methods utilize the Weibull distribution form or the value of β .

· Method 1 (the sample median) does not assume any distribution form.

· Method 2 (MLE, 2 variables) assumes the Weibull distribution form of Z but leaves shape parameter β as an unknown variable.

· Method 3 (MLE, 1 variable) both assumes the Weibull distribution form and uses the value of shape parameter β .

· Method 4 (adaptive Bayesian) both assumes the Weibull distribution form and uses the value of shape parameter β .

The situation where the actual distribution deviates from the assumed model distribution will be investigated in Section 5.

The relative errors for β = 2.5 are reported in

Accordingly, MLE (1 variable) has the smallest relative error since it assumes both the correct distribution form and the correct shape parameter. Method 4 is based on binary outcomes (escape or no escape) observed in tests using the method of constant stimuli. Intuitively, a binary outcome provides less information than a full sample of escape time. The relative error of method 4 is larger than that of method 3 even though they have the same assumptions. However, we will learn in Section 5 that method 4 is very robust with respect to an incorrect model while for method 3, an incorrect model is fatal. The relative error of method 4

Method | Relative error of inferred z ( m ) | |||
---|---|---|---|---|

N = 2 | N = 20 | N = 40 | N = 80 | |

Sample median | 31.21% | 12.43% | 8.95% | 6.39% |

MLE (2 variables) | 32.26% | 10.47% | 7.43% | 5.25% |

MLE (1 variable) | 28.30% | 8.94% | 6.34% | 4.47% |

Adaptive Bayesian, prior σ z / μ z = 40 % | 29.00% | 12.45% | 9.00% | 6.42% |

Adaptive Bayesian, prior σ z / μ z = 20 % | 17.94% | 10.83% | 8.31% | 6.15% |

Adaptive Bayesian, prior σ z / μ z = 10 % | 9.71% | 7.88% | 6.72% | 5.41% |

is comparable to that of method 1. This feature will be examined in more detail in Section 5. An interesting result of

The relative errors for β = 5 are reported in

We study numerically the relative error ( err ( k ) ) of each inference method as the number of trials (N) is increased. The simulation results for β = 2.5 are shown in

Method 4 (adaptive Bayesian) updates the prior distribution based on binary

Method | Relative error of inferred z ( m ) | |||
---|---|---|---|---|

N = 2 | N = 20 | N = 40 | N = 80 | |

Sample median | 16.03% | 6.23% | 4.48% | 3.20% |

MLE (2 variables) | 16.28% | 5.24% | 3.72% | 2.62% |

MLE (1 variable) | 15.17% | 4.50% | 3.18% | 2.24% |

Adaptive Bayesian, prior σ z / μ z = 40 % | 22.05% | 6.71% | 4.67% | 3.27% |

Adaptive Bayesian, prior σ z / μ z = 20 % | 14.69% | 6.30% | 4.53% | 3.22% |

Adaptive Bayesian, prior σ z / μ z = 10 % | 8.98% | 5.45% | 4.17% | 3.08% |

outcomes observed in the method of constant stimuli. Unlike the three methods in the left panel, method 4 has a finite relative error (the prior) before any binary outcome is observed. Intuitively, the prior can be viewed as an existing set of n 0 binary outcomes. A narrower prior corresponds to a larger existing set (larger n 0 ). Consequently, the relative error is expected to be proportional to 1 / N + n 0 .

err ( N ) = C 4 N + n 0 , n 0 = ( C 4 err ( prior ) ) 2 , err ( prior ) = σ z μ z (29)

The right panel of

We explore how the relative error changes as the shape parameter β of the Weibull psychometric function is varied. We reiterate that in all simulations of Section 4, both the actual and the model psychometric function have the same Weibull form and the same shape parameter β . The situation where the actual psychometric function deviates from the model will be studied in the next section.

As demonstrated in

Bayesian method. We examine the behaviors of relative errors in several situations where the model deviates from the actual psychometric function.

When both the actual psychometric function and the model have the Weibull form, we write them as two functions of x = log 10 u with different shape parameters and medians.

F actual ( x ) = ϕ ( x − y ( m , e ) ; β ( e ) )

F model ( x ) = ϕ ( x − y ( m ) ; β )

where function ϕ ( s ; β ) is defined in (18). In the actual psychometric function, β ( e ) is the true shape parameter and y ( m , e ) the true median of log 10 Z . In the model psychometric function, β is the prescribed shape parameter and y ( m ) the inference variable.

The results for the 4 cases of ( β , β ( e ) ) = { ( 2.5,2.5 ) ; ( 2.5,5 ) ; ( 5,2.5 ) ; ( 5,5 ) } are displayed in

· the β value cannot be set arbitrarily in the MLE (1 variable); and

· when the true β ( e ) of the actual psychometric function is unknown, the method of MLE (1 variable) is simply not applicable.

One striking feature of the adaptive Bayesian method is that it converges to the true median even when a wrong β value is prescribed in the model (right panel of

When the prescribed β in F model ( x ) is larger than the true β ( e ) in

F actual ( x ) , the model distribution is narrower than the actual distribution. When β deviates from β ( e ) in the direction of β > β ( e ) , the relative error becomes noticeably larger than that of β = β ( e ) . In addition, the decay of relative error vs N no longer follows the trend of 1 / N + n 0 . Fitting C / ( N + n 0 ) d to the relative error vs N for { β = 5 , β ( e ) = 2.5 } yields an exponent of d = 0.43 instead of 0.5. The fitting curves respectively for β > β ( e ) and for β = β ( e ) have the expressions

err = 0.50 ( N + 1.6 ) 0.43 for β = 5 > β ( e ) = 2.5

err = 0.58 N + 2.1 for β = β ( e ) = 2.5

When the prescribed β in F model ( x ) is smaller than the true β ( e ) in F actual ( x ) , the model distribution is wider than the actual distribution. When β deviates from β ( e ) in the direction of β < β ( e ) , the relative error increases only slightly from that of β = β ( e ) . The relative error vs N is still approximately described by the law of 1 / N + n 0 , with the coefficient slightly above that of β = β ( e ) . The fitting curves respectively for β < β ( e ) and for β = β ( e ) are

err = 0.35 N + 0.74 for β = 2.5 < β ( e ) = 5

err = 0.30 N + 0.56 for β = β ( e ) = 5

In the right panel of

We consider the situation where the model distribution of subjective threshold Z is Weibull while the actual distribution of Z is gamma, which has the form

ρ ( gamma ) ( z ) = 1 Γ ( κ ) θ ( z θ ) κ − 1 exp ( − z θ ) (30)

where κ is the shape parameter and θ the scale parameter of the gamma distribution. Two Weibull distributions for β = { 2.5 ; 5 } and three gamma distributions for κ = { 4 ; 8 ; 16 } are compared in the left panel of

F actual ( x ) = 1 Γ ( κ ) γ ( κ , 10 x θ )

F model ( x ) = ϕ ( x − y ( m ) ; β )

where γ ( κ , s ) is the lower incomplete gamma function defined as

γ ( κ , s ) = ∫ 0 s t κ − 1 e − t d t

As in the case of Weibull distribution, when the actual distribution of Z is gamma, the relative error of inferred z ( m ) ≡ med ( Z ) is independent of the scale parameter θ .

We test the adaptive Bayesian method using the Weibull model with β = 2.5 on three actual psychometric functions of gamma distributions with respectively κ = 4 , 8, and 16. Relative errors vs N are plotted in the right panel of

Now we look at the most important reason why we choose to infer med ( Z ) instead of E ( Z ) . In the assumed model (Weibull), the mean is related to the median by

E ( Z ) = Γ ( 1 + 1 β ) ( ln 2 ) − 1 β med ( Z ) (31)

which varies with β . In the actual psychometric function, the relation between

the mean and the median is in general different from (31). Given the convergence of inferred med ( Z ) , it follows that if we choose to infer E ( Z ) , then the inferred E ( Z ) does not converge to the true E ( Z ) . That is, if we choose to infer E ( Z ) , the convergence of inference is lost.

In the left panel of

The results indicate that it is more important to accommodate the wider end in a possible range of the actual distribution width. If the actual distribution has a large width, failing to accommodate it will lead to a significant increase in relative error (compare triangles and diamonds in

In the adaptive Bayesian method, the posterior distribution of y ( m ) = log 10 z ( m ) not only produces a posterior median as the inference result for y ( m , e ) , but also

yields a posterior standard deviation σ y ( post ) . When the sample size N is not very small, the posterior of y ( m ) is approximately a normal distribution. A normal posterior for y ( m ) translates to a log-normal posterior for z ( m ) . We view the posterior standard deviation of z ( m ) divided by the posterior median of z ( m ) as the predicted relative error in the inferred z ( m ) . This prediction is practically operational since it comes out naturally from the posterior. The predicted relative error is expressed in terms of σ y ( post ) using (26).

( Predicted relative error in inferred z ( m ) ) ≡ σ z ( post ) μ z ( post ) = exp ( q 0 2 ( σ y ( post ) ) 2 ) [ exp ( q 0 2 ( σ y ( post ) ) 2 ) − 1 ] , q 0 = ln 10 (32)

The blue line and red line in the left panel of

Thus, using a wide model distribution yields a predicted error that is too conservative. A wide model distribution is less a problem for minimizing the actual error than for precisely predicting the error. Synthesizing the results of subsections 5.2 and 5.3, we draw two conclusions regarding the adaptive Bayesian method:

· when presented with the choice of two likely candidate model distributions, to minimize the actual error of inference in this uncertain situation, the best strategy is to select the wider one; and

· the error predicted from the posterior standard deviation is dominated by the assumed model distribution and is unreliable; the actual error of inference is mainly influenced by the actual distribution of the subjective threshold.

The sample median is calculated from independent samples of the subjective threshold Z measured in the method of limits. The adaptive Bayesian method utilizes binary outcomes observed in the method of constant stimuli. Intuitively each binary outcome contains less information than a full sample value of Z: it only indicates whether the value of Z in that trial is below or above the prescribed stimulus. Given the difference in information content between these two types of data, we like to find out if the adaptive Bayesian method based on binary outcomes is less efficient in inferring med ( Z ) than the sample median calculated from full samples of Z.

from the actual psychometric function, the Bayesian method is as efficient as the sample median in estimating med ( Z ) . For moderately large N, the relative errors of the Bayesian method are indistinguishable from those of the sample median, and are fairly independent of the prior. Thus, we conclude that for estimating the point of subjective equality (PSE) in psychophysical experiments, the method of constant stimuli is as efficient as the method of limits.

We considered the psychophysical experiments in which a skin area of a test subject is exposed to an electromagnetic heating or a contact heating. Once the heat induced exceeds the subject’s personal pain tolerance, the subject escapes from the heating. In this stimulus-response problem, when the specifications of heating (the spot size and power density of the electromagnetic beam) are fixed, the stimulus is solely described by the exposure duration, and the escape response is governed by a subjective threshold that varies among different subjects and varies among different realizations for the same subject. Mathematically, the subjective threshold is a random variable (Z). In the context of psychometric function, the median subjective threshold is the point of subjective equality (PSE). We studied the problem of inferring the median subjective threshold in psychophysical experiments. There are two types of experiments: a) the method of limits, and b) the method of constant stimuli. In the method of limits, the heating is kept on until escape (i.e., the stimulus keeps increasing until escape). The recorded escape time gives a sample value of the subjective threshold for that particular trial. In this way, independent samples of the subjective threshold are measured. In the method of constant stimuli, the heating is kept on only for a prescribed duration (i.e., the stimulus is prescribed). The observed outcome is binary: escape or no escape. In the inference, we modeled the subjective threshold as a Weibull distribution which is described by shape parameter β and scale parameter α . The actual distribution of the subjective threshold may deviate from the assumed model.

We examined four methods for inferring the median subjective threshold: 1) sample median, 2) maximum likelihood estimation (MLE) with 2 unknown variables ( α and β are both unknown), 3) MLE with 1 unknown variable ( β is given and α is unknown), and 4) adaptive Bayesian method. Methods 1 - 3 are based on independent samples of subjective threshold measured in the method of limits. Out of methods 1 - 3, method 1 does not assume any distribution form; not surprisingly it has the largest inference error. Method 3 assumes both the Weibull distribution form and the exact value of β ; it has the smallest inference error. However, when the assumed β value is wrong or the actual distribution deviates from the assumed Weibull model, method 3 produces an incorrect result.

A large part of our study was focused on the adaptive Bayesian method, which assumes a Weibull model distribution and uses it to extract information from binary outcomes observed in the method of constant stimuli. Based on the Monte Carlo simulation results, we draw several key conclusions about the adaptive Bayesian method.

· Even when the actual psychometric function deviates from the assumed model, the inferred med ( Z ) in the adaptive Bayesian method converges to the true value.

· If we choose to infer E ( Z ) (the mean subjective threshold) in the adaptive Bayesian method, this nice convergence property is lost.

· At any fixed N (number of trials), the actual relative error of the adaptive Bayesian method is mainly influenced by the width of the actual distribution of Z: a wider actual distribution leads to a larger relative error.

· When the width of the assumed model distribution is wider than that of the actual distribution, the actual relative error is fairly insensitive to the inconsistency. In contrast, when the width of the assumed model distribution is narrower than that of the actual distribution, the actual relative error increases significantly with the inconsistency.

· The predicted relative error from the posterior standard deviation is mainly influenced by the width of the assumed model distribution. Using a model distribution that is too narrow produces both a large actual error and a false sense of accuracy (small predicted error). On the other hand, using a model distribution that is too wide yields a predicted error that is too conservative.

· When facing an uncertainty in the actual distribution width, the best strategy is to select a smaller β in the Weibull model to accommodate the wider end in a possible range of the actual distribution width.

· The actual errors of the adaptive Bayesian method are indistinguishable from those of the sample median, and are fairly independent of the prior used.

In summary, the adaptive Bayesian method is as efficient as the sample median for estimating the point of subjective equality in psychophysical experiments. This efficiency of the Bayesian method is remarkable since the sample median requires measuring full sample values of the subject threshold in the method of limits. In contrast, the adaptive Bayesian utilizes the binary outcomes observed in the method of constant stimuli. In addition, the Bayesian method achieves this efficiency even when the actual psychometric function deviates from the assumed model. In the adaptive Bayesian method, the error predicted from the posterior standard deviation is unreliable; it does not reflect the actual error. The actual error is mainly determined by the actual distribution of Z while the predicted error is dominated by the assumed model. To minimize the actual inference error when the actual distribution width is uncertain, it is best to use a wider model distribution.

The authors acknowledge the Joint Intermediate Force Capabilities Office of U.S. Department of Defense and the Naval Postgraduate School for supporting this work. The views expressed in this document are those of the authors and do not reflect the official policy or position of the Department of Defense or the U.S. Government.

The authors declare no conflicts of interest regarding the publication of this paper.

Wang, H.Y., Adamzadeh, M., Burgei, W.A., Foley, S.E. and Zhou, H. (2021) Inference of Median Subjective Threshold in Psychophysical Experiments. Journal of Applied Mathematics and Physics, 9, 982-1002. https://doi.org/10.4236/jamp.2021.95068