A novel paradigm to test regression

We designed a novel psychophysical experiment in which participants had to extrapolate a noisy parabola displayed on a computer screen. In each trial, we chose the parameter w of the parabola y = wx² from a bimodal prior distribution π(w) where the two modes are centered at w = 1 and w = −1 and the variances are given by (see Eq 5). The parameter w was either positive (parabola facing upwards) or negative (parabola facing downwards), with the same probability, i.e., 0.5. We selected four dots on the parabola with x-positions close to the parabola’s vertex and added zero-mean Gaussian generative noise σ_g to the dots’ y-positions (see Eq 4). We then presented a fifth dot to the right of the stimulus, always at the same x-position x^⋆ = 2. Participants could move the fifth dot up and down along the y-axis by using the up and down arrow keys. Participants were asked to adjust the y-position so that the dot correctly extrapolated the parabola. During the adjustment task, participants saw only the 4-dot stimulus but not the generating parabola. After the the participant had validated his/her response, we showed the generating parabola and the adjusted point as feedback. Participants were naive to the purpose of the study. They were not informed about the existence of a prior distribution of the parabola’s quadratic parameter, the parabola’s bimodality nor the level of generative noise.

In our main experiment, we set the standard deviation of each prior mode to σ_π = 0.1 (if not specified otherwise, assume this value throughout this work) and fixed the values of x-positions to x₁ = −0.3, x₂ = −0.1, x₃ = 0.1 and x₄ = 0.3. We generated j ∈ (1, …, 20) unique stimuli at a low (0.03), medium (0.1) and high (0.4) value of the generative noise σ_g. The rationale is that the higher the noise level, the higher the uncertainty (the lower the likelihood) and the more participants rely on the prior if they act consistently with a Bayesian regression model. At each noise level, we ran 400 trials, repeating each unique stimulus D_j 20 times. The stimulus presentation order was randomized within each noise level (Fig 1B). Thus, we obtained a set of responses for each noise level and for each of the 20 stimuli . The advantage of observing several responses for the same exact stimulus is that we can compare the observed response distribution to the response distributions predicted by the different models.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Experimental protocol.

(A): Procedure of a single trial. First, a fixation dot was presented for 1s before the 4-dot stimulus appeared. Observers then had unlimited time to adjust the fifth dot with the up and down arrow keys. They then clicked the space bar to confirm the final position of the adjustable dot. After the response, the generative parabola was shown for 1s as feedback. (B): Experiment 1: The experiment consisted of two sessions on two separate days. Both sessions began with 10 practice trials with virtually no noise (σ_g = 10⁻⁵), followed by 4 blocks of 50 trials of low noise (σ_g = 0.03). In session 1, the low noise blocks were followed by 8 blocks of 50 trials of medium noise (σ_g = 0.1), while in session 2, the low noise blocks were followed by 8 blocks of 50 trials of high noise (σ_g = 0.4). In total, each participant completed 400 trials per noise level, with 20 repetitions of 20 unique stimuli. In this experiment, σ_π was set to 0.1. (C): Experiment 2: the experiment consisted of a single session which began with 20 practice trials with very low noise (σ_g = 10⁻²), followed by 10 blocks of medium noise (σ_g = 0.1) trials. Each block consisted of 20 trials, with the generative parabola shown as feedback, as in Experiment 1. Half of the 200 trials consisted of stimuli which were presented just once, while the remaining 100 trials consisted of 10 repetitions of 10 unique stimuli. In this experiment, σ_π was set to 0.5. See Materials and methods for more details.

https://doi.org/10.1371/journal.pcbi.1007886.g001

Seven naive participants took part in the experiment. We denote the set of all responses of participant k by . The stimulus presentation order was identical for each participant. At the beginning of each experimental session, we showed a virtually noiseless 4-dot stimulus σ_g = 10⁻⁵ to familiarise the participants with the task and to estimate their internal noise sources (explained in more detail below).

In a second experiment, we set the prior standard deviation to σ_π = 0.5 and the generative noise to the medium level of σ_g = 0.1. Instead of using fixed x-positions for the stimuli we added weak Gaussian noise (see Material and methods). We repeated 10 unique stimuli 10 times each, which yielded a total of 100 trials. Four naive participants (different from those recruited for the first experiment) completed the experiment. The key difference to the σ_g = 0.1 condition in the first experiment was that here, the prior provided much less information about the curvature of the parabolas.

In total, we studied four different conditions: three with σ_π = 0.1 and σ_g ∈ {0.03, 0.1, 0.4} in experiment 1, and one with σ_π = 0.5 and σ_g = 0.1 in experiment 2.

The regression models

We considered five regression models (see Materials and methods). The Maximum likelihood regression (ML-R) model computes only the point estimate of w that maximizes the likelihood p(D_j|w) and does not make use of the prior distribution at all. The maximum a posteriori model (MAP-R) combines the likelihood with the prior distribution to compute the mode of the posterior distribution p(w|D_j). Despite the fact that it uses the bimodal prior, MAP-R cannot produce bimodal responses because it relies on a point estimate of w. The Bayesian regression (B-R) model takes the entire posterior distribution into account: (2) In B-R, the noise level σ_g plays the crucial role of modulating the relative strength of unimodal likelihood and bimodal prior, and hence determines the transition between a unimodal and a bimodal response distribution. Relaxing the assumption that the true generative noise is known, we included an additional variant of B-R that (deterministically) estimates the generative noise for the current stimulus D_j. We indicate this variant of B-R with a subscript: B-R_σ. Note that technically, this trial-by-trial noise estimation could also be applied to the other models such as MAP-R. However, MAP-R reduces the posterior distribution over the quadratic parameter to a point estimate. Therefore, an additional trial-by-trial noise estimation would not change its prediction substantially, i.e., it would only shift the unimodal prediction but would not induce a bimodal predictive distribution. Thus, we did not consider a corresponding “MAP-R_σ” variant. As a null model, we included prior regression (P-R), which replaces the posterior with the prior, i.e., it does not use the likelihood.

For all models considered here the predictive distribution depends deterministically on the 4-dot stimulus. In this sense, they rely on exact inference. Noisy inference is an alternative which assumes that the inference process is corrupted by noise [29]. This alternative would require an additional noise parameter which governs the level of inference noise (see S1 Text). Here, we constrain ourselves to models with exact inference to remain fitting-free, i.e., model predictions for a given stimulus D_j have no free parameter. We used the true values of the hyperparameters because we assume that the participants learned the generative model within a few trials (see S1 Text). Thus, the model predictions require no fitting. For more details, see Materials and methods.

In the plots, we denote the models by the arguments of their predictive distributions, i.e., y|x, w_ML for Maximum Likelihood regression (ML-R); y|x, w_MAP for Maximum a Posteriori regression (MAP-R); y|x, D for Bayesian regression (B-R); y|x, D, σ_g for Bayesian regression with noise estimation (B-R_σ) and y|x for prior regression (P-R).

The decision models

We considered two decision models that turn the predictive distributions into a response distribution: probability matching and Bayesian decision theory. In the case of probability matching, it is assumed that participants draw random samples from the predictive distribution: . If not stated otherwise, we use sampling-based decisions throughout this work. Meanwhile, according to Bayesian decision theory, participants select a response by minimizing the expected loss function . Here, we considered only the square loss, which is equivalent to the choosing mean of the predictive distribution . Independently of the form of the loss function, the Bayesian decision theory generates responses from the predictive distribution deterministically. When we use loss-based decision models, we indicate this by adding the prefix “L”: to the model, e.g., L: y|x, D for a loss-based decision model applied to Bayesian regression.

To model participants’ responses, we also accounted for internal sources of noise, i.e., noise which is inhere to neural processing, decision making and the execution of motor action [29, 39]. We call the sum of these noise components motor noise for brevity. The motor noise is not a model parameter but a participant-specific parameter. We computed the motor noise σ_m for each participant from the 20 responses to the noiseless stimulus. To ensure robustness to outliers, we used the average value between the 16% and 84%-percentile of the response distribution. The values of motor noise for the seven participants of the main experiment were , respectively while for the second experiment, we used the average of these values, i.e., 0.48 because the noise-free responses of experiment 2 were not available. The motor noise was included in the models by convolving the predictive distribution with a Gaussian of variance . In the case of loss-based decision models, motor noise was the only source of response variability.

Modality of predicted and observed response distributions

Fig 2(A) shows the responses of a representative participant along with the predicted response distributions of the different models. Both ML-R and MAP-R ignore one of the modes (here, the mode corresponding to a downward-facing parabola). In addition, the parabola predicted by ML-R has lower curvature than the parabolas predicted by any of the other models (i.e., the absolute value of the ML-R parabola’s quadratic parameter is lower) than that of the parabola that the participant responded with. A potential explanation for this finding is that, while ML-R does not take the prior into consideration, humans do make use of the prior. In the low noise regime (σ_g = 0.03, Fig 2(B)), the discrepancy between the participant’s response distribution and the prior regression model’s predicted response distribution in terms of the number of modes (unimodal and bimodal, respectively) rules out the explanatory validity of the latter. In the higher noise regimes (σ_g ∈ (0.1, 0.4), Fig 2(C) and 2(D)), MAP-R and ML-R fail to account for the fact that the participant’s responses are distributed across both modes. The finding that at σ_g = 0.03 MAP-R matches the participant’s responses often with high accuracy provided implicit evidence that participants used the prior and had learned the parabola’s generative model.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Example responses.

B-R is the only model that can explain the transition from unimodal response (at low noise, (B)) to bimodal response distribution (at high noise, (D)). (A) A sample stimulus (green dots) at high noise level (σ_g = 0.4). For this specific stimulus, contours indicate the response distributions predicted by ML-R, MAP-R, B-R and B-R_σ (not shown to the participant) at various x^⋆. At x^⋆ = 2, we recorded the participant’s responses (gray dots). The cross section at x^⋆ = 2 is shown in (D). (B—E) The predicted response distributions at x^⋆ = 2 of ML-R (blue), MAP-R (orange), B-R (red), B-R_σ (dark red), P-R (green) and observed responses (gray). As σ_g increases (B—D), the data becomes less informative. Consequently, and in accordance with B-R, the response distribution becomes more bimodal. (E) Due to the weak prior the predictions of B-R and B-R_σ respond more strongly to the data and diverge from the modes of P-R more stronlgy than in the previous conditions. The skewness of B-R_σ results from the mixture of both Gaussian components.

https://doi.org/10.1371/journal.pcbi.1007886.g002

Fig 2(E) illustrates the participant’s responses in the condition when σ_π = 0.5. The participant’s responses cover a wider range of values than in the conditions of experiment 1 when σ_π is smaller (i.e. σ_π = 0.1). While the generative noise σ_g = 0.1 is the same as in Fig 2(B), this condition is more difficult because the prior is less reliable. As a consequence participants rely more strongly on the noisy stimulus and produce more response variability. In this example, the responses are closely clustered around the center. B-R_σ is attracted more strongly to the center than B-R because the former is more driven by the stimulus due to underestimating the noise.

B-R outperforms the other models

Fig 2. In order to formally assess model performance, we next conducted a quantitative model comparison across all participants. For each of the seven participants individually, we computed the log probability that the participant’s responses arise from the given model. We summed these log probabilities for all of the unique stimuli D_j as a measure of the quality of the model. Fig 3(A) shows these values relative to the B-R baseline value for each noise level, averaged across participants. Negative values indicate poor performance relative to B-R. A subject-level analysis showed that the model comparison results were not driven by any single participant’s data (see S1 Text). Because our model comparison is fitting-free, we do not need to account for different levels of model complexity. Indeed, in the present case, the log likelihood comparison is equivalent to using the Bayesian Information Criterion.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Model comparison.

The model comparison shows that the B-R model best explains the data (A, B) and that sampling-based decision models outperform loss-based decision models (C). (A) Difference in log likelihood with respect to B-R averaged over participants for different experimental conditions. Negative values mean that B-R wins the comparison. B-R is either winning (σ_g ∈ (0.03, 0.1)) or equivalent to P-R because the two coincide at high levels of parameter uncertainty (σ_g = 0.4 and σ_π = 0.5). (B) The expected likelihood of each model for a randomly selected participant shows what fraction of participants are best described by a model. Overall, B-R and B-R_σ describe the population best. (C) Log likelihood difference between a sampling and a loss-based decision model. Negative values favour sampling. At all other conditions and for all regression models, sampling explains the data better than loss-based decision models with exact inference. For B-R, B-R_σ and P-R, loss-based models do not predict bimodal responses. At low noise σ_g = 0.03, loss-based models underestimate the response variance. Error bars represent the SEM across participants.

https://doi.org/10.1371/journal.pcbi.1007886.g003

As Fig 3(A) shows, B-R is among the highest performing models for all conditions. As the task difficulty increases (left to right), P-R performance approaches that of B-R. This is because the parameter uncertainty encoded in the prior becomes more important and the response distribution becomes bimodal. Neither MAP-R nor ML-R can capture this and therefore perform poorly. These results are consistent across participants (see S1 Text for a subject-level analysis).

At low noise σ_g = 0.03, participants give unimodal answers and the mean predictions of B-R and MAP-R are indistinguishable. Then the model that better captures the response variability wins. In general, B-R explains the variability of the responses better. The variability is also the reason why P-R performs relatively well under the more difficult conditions, i.e., σ_g = 0.03 and σ_π = 0.5.

The averaged results are largely consistent with a subject-level analysis. A notable exception is that at σ_g = 0.03, MAP-R emerges as the best model (closely followed by B-R) for participants 3, 4 and 7 (see S1 Text). A Bayesian random effects analysis confirms this. Specifically, we used the model posterior averaged over a randomly selected participant k. This measure reflects the ratio of participants for which model wins. Fig 3(B) shows that the responses of the majority of participants are best modelled by B-R or B-R_σ. Since B-R interpolates between MAP-R (at low noise) and P-R (at high noise), as expected, at σ_g = 0.03, i.e., the easiest condition, the responses of some participants are also well modelled by MAP-Rwhile at the most difficult condition, i.e., when σ_π = 0.5, the responses of half of the participants are best described by P-R (and the other half by B-R).

Next, we investigated if sampling or the loss function perspective explains the responses better. Fig 3(C) depicts the log likelihood of loss-based decision making compared to sampling for each model. Negative values indicate that sampling wins. Sampling explains the data better for all models and in all experimental conditions. One explanation is that the loss mechanism turns bimodal predictive distributions into unimodal predictive distributions. Here, we use the square loss such that the (unimodal) response distribution is centered on the mean of the predictive distribution. In the case of P-R, the mean of the predictive distribution lies at the center of both modes. Clearly, this method does not capture bimodal responses. This is why the performance difference between the two decision models is smallest at σ_g = 0.03, where all models except for P-R make predictions which are close to unimodal.

The second explanation for the better performance of sampling is that the loss function approach with exact inference underestimates response variability. The response variability differs from one stimulus to another and is often higher than σ_m. This explains the better performance of sampling for MAP-R and ML-R, since the effect of turning a bimodal response distribution into a unimodal one is absent. In these cases, the sampling-based decision model has the effect of increasing the variance of the predicted response distribution by . This leads to better model performance on variable response data, even in the experimental conditions σ_g = 0.03 where participants respond unimodally.

In conclusion, from the models considered here, B-R with sampling best explains participants’ responses.

B-R explains the generative noise-dependent increase in response variance

A key characteristic of B-R is the transition of the model’s posterior predictive distribution from unimodality to bimodality as σ_g increases, i.e., as the data become less informative. To analyse this transition, we used the participants’ response variances at the different levels of generative noise.

The variance of the response distributions is sensitive to bimodality. For example, if all responses are distributed evenly across both modes, the variance is close to 16, which corresponds to the variance of the prior P-R. If all responses are located in a single mode, the variance is typically smaller by a factor of ten (e.g., see Fig 2(B)). We explain this in more detail below.

For each stimulus D_j and for each participant, we computed the variance of the 20 observed responses . Because we have 7 participants, this yields a distribution over 7 × 20 = 140 empirical variance values at each value of σ_g. We compared this distribution with the response variance distribution predicted by the models. To achieve a higher resolution and show the dynamics of the variance as a function of the generative noise, we generated 5000 unique stimuli from a densely spaced σ_g instead of relying on the small number of stimuli and noise levels used in the experiment (see Materials and methods).

The empirical variance distribution (gray) and its median (black) are shown in Fig 4 along with the predicted median of the variance distribution for each model (color coded). For B-R and B-R_σ, we plotted the distribution in Fig 4(A) and 4(B), respectively.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Response variances of predicted and empirical distributions, as a function of generative noise.

B-R best explains the increase in response variance as a function of the generative noise σ_g. Variances of the empirical response distributions from all participants (gray dots, median: gray line) and predicted response distributions, corresponding to the two B-R variants (median: red line, log probability: heatmap). B-R (A) Interpolating between MAP-R and P-R, only the B-R variants capture the upward trend in the data. At σ_g = 0.4, B-R fails to account for the empirical responses with close-to-zero variances. (B) At σ_g = 0.4, B-R_σ predicts a bimodal variance distribution because, in trials with low noise estimates, the predicted response distribution is unimodal and thus variance is low. Because of these low-variance trials, the median of B-R_σ increases slower than the median of B-R and captures the empirical median better. Because ML-R and MAP-R behaved identically, the MAP-R represents both regression models.

https://doi.org/10.1371/journal.pcbi.1007886.g004

The median of the B-R variance distribution increases with the noise level. This is due to the fact that B-R’s predicted response distribution transitions from unimodal to bimodal; this transition is modulated by the generative noise, which determines the relative contribution of the prior and likelihood to the response distribution. Consequently, the B-R model is the only one for which the variance values smoothly transition from the MAP variance at the low noise level (σ_g = 0.03) to the P-R variance at the high noise level (σ_g = 0.4). The P-R variance remains constant and the MAP-R variance increases very weakly as a function of the generative noise. The variance analysis provides further evidence for the superiority of both B-R variants over the other models.

While B-R captures the general trend in the data, it fails to account for two key characteristics. First, the median variance increases slower than B-R would predict, and secondly, at high levels of generative noise, B-R fails to reproduce the lower part of the distribution (where response variances are close to zero). A potential explanation for this discrepancy is that participants estimate the noise on a trial-by-trial basis. When the noise added to the 4-dot stimulus was, by chance, such that the dots appeared to be well-aligned on a parabola, participants would, presumably, underestimate the generative noise and respond in a way which was consistent with a unimodal distribution. The fact that the B-R_σ model captures the empirical variance better than B-R provides some evidence for this idea. On some trials, B-R_σ underestimates the true σ_g and applies the B-R formalism with high confidence in the stimulus data. In these cases, the model relies strongly on the likelihood and bimodality, which normally enters through the prior, is not achieved. Rather, the resulting response distribution is unimodal and has low variance.

Despite the fact that B-R_σ describes the qualitative features of the variance distribution better than B-R, it performs worse in terms of log likelihood. This shows that the low variance responses of humans and of B-R_σ do not always coincide on a trial-by-trial basis.

To better understand the relation between response variance and bimodality, we dissect the variance of a bimodal response distribution to stimulus D_j into its components: (3) where μ₁ and μ₂ are the means of the modes of the posterior predictive distribution, c is the mixture coefficient and corresponds to the variance of both modes. The unimodal contribution is not mode-specific because we chose a symmetrical prior. The first two terms constitute a unimodal contribution and the last term a bimodal contribution. The latter is controlled by the mean dispersion (μ₁ − μ₂)² and a prefactor c(1 − c) that is equal to zero for c ∈ {0, 1} and is maximal for c = 1/2. To determine to what extent each component of this dissection is present in the response data, we defined the empirical counterparts of μ₁, μ₂ as the means of the upper and lower modes of the response distribution and c as the mixture coefficient, corresponding to the fraction of positive responses r > 0. For the unimodal variance contribution , we used the variance of the mode which contains the majority of responses (see Methods for more details). The comparison between data and models shows that both B-R variants correctly predict the driver of the observed variance to be the transition to bimodality. Fig 5(A) shows the predicted positive coefficient c (median) of the models as a function of the empirically-observed coefficient across all participants and stimuli (in the main experiment). As further evidence for the validity of the B-R perspective, both B-R variants correctly predict the fraction of positive responses. Indeed, the smooth transition from a unimodal to a bimodal distribution is nicely captured by B-R. In contrast, MAP-R transitions sharply and is more reminiscent of a step function while P-R predicts equally strong modes across all noise level conditions.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Median of each of the bimodal response distribution variance components across all participants and stimuli.

(A) Predicted coefficient of positive mode as a function of the empirical coefficient (across all noise levels). ML-R behaves identically to MAP-R. Thus, the MAP-R curve represents both models. The shaded area shows the 40% and 60% quantiles. (B) Prefactor of bimodal contribution as a function of generative noise. Data jittered for visibility. (C) Unimodal contribution to the variance. Empirical variance computed on mode with majority of responses. (D) Mean dispersion. Only trials with bimodal responses included. As the stimulus becomes more noisy, human responses and B-R variants conform to the prior.

https://doi.org/10.1371/journal.pcbi.1007886.g005

The bimodal distribution of responses depends on the prefactor c(1 − c) and the mean dispersion. Fig 5(B) shows the median value of the prefactor as a function of generative noise (across all participants and stimuli). Data and model predictions qualitatively match the behaviour of the variance in Fig 4. Indeed, the other contributions to the variance are less important. The mean dispersion, shown in Fig 4(C), plays the role of a large constant. The unimodal contribution to the variance, shown in Fig 5(D), is small compared to the bimodal contribution. In conclusion, the coefficient c plays the dominant role in determining the variance of the response distribution. Because the two B-R variants estimate c sufficiently well, they best match the empirical variance distribution.

Interestingly, all models overestimate the unimodal variance, with the exception of MAP-R in the low noise condition (Fig 5C). The B-R variants predict larger variance than MAP-R because they translate the posterior parameter into response uncertainty. P-R predicts even larger response variance because it uses the prior parameter uncertainty which is generally larger than the posterior one. Despite the fact that MAP-R best describes the median variance, it performs worse than B-R in terms of log likelihood. Fig 5(C) reveals that one factor contributing to the poorer performance of MAP-R is the occurrence of unimodal, high variance responses.

In summary, the variance analysis provides further evidence that B-R captures the way in which generative noise induces a transition from unimodality to bimodality in participant responses. However, B-R overestimates response variance. Trial-by-trial estimation of the noise offers a potential explanation for why participants cluster their responses more unimodally than predicted.

Unimodal responses are overall best explained by B-R

Thus far, the main factor behind the superiority of the B-R model’s performance relative to the other models is the ability of B-R to capture the bimodality of responses, i.e., to correctly set the mixture coefficient. However, the previous analysis showed that unimodal variance decreases as a function of generative noise while B-R predicts an increase. It remains unclear if B-R still wins the model comparison in a unimodal setting where performance is independent of the mixture coefficient.

To address this question, we conducted a model comparison on a unimodally conditioned dataset. For each stimulus D_j, we considered only responses in which the dominant response mode and the dominant mode of the model predictions coincided (see Methods for details). The conditioning yields a unimodal dataset in the sense that all predictions and responses belong to the same mode. To make the model comparison fair for the bimodal predictive distributions of P-R and the B-R variants, we removed the inferior mode and normalised the remaining probability mass to one.

At the group level, B-R wins the model comparison across all conditions, as shown in Fig 6(A). B-R clearly outperforms the other models in two conditions in particular: at σ_g = 0.03 and at σ_π = 0.5.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Unimodal model comparison.

The unimodal analysis confirms previous results: overall B-R with sampling wins the model comparison. (A) Differences in log likelihood on unimodal data, averaged over participants. Negative values mean that B-R wins. ML-R is omitted because its poor performance complicates visualisation. (A) B-R wins at σ_g = 0.03 and σ_π = 0.5, but not in the other conditions. (B) All models use the quadratic loss function to select responses, with response variance given by the motor noise . B-R with sampling explains the unimodal data best for most participants. High subject-level variability results in large errors (see S1 Text for a subject-level analysis). (C, D) The fraction of participants best described by a given model. At σ_g = 0.4, several models perform well. Error bars indicate SEM across participants.

https://doi.org/10.1371/journal.pcbi.1007886.g006

Interestingly, in the bimodal dataset, B-R did not emerge as a clear winner at σ_π = 0.5 because P-R performed similarly well. Thus, B-R is not better than P-R at modelling how participants balance the two modes, but once the mode is chosen, it performs better. At σ_g = 0.1, B-R and B-R_σ outperform the other models. B-R_σ wins by a small margin (see axis scaling). However, a subject-level analysis (see S1 Text) shows that the average is mostly driven by participant 1, while in the case of other participants all models perform similarly well. At σ_g = 0.4, no clear winner emerges. Intuitively, this makes sense because the stimulus is not informative and all models rely mostly on the prior information about w. Here, B-R wins or performs similarly to other models.

The Bayesian random effects analysis (results shown Fig 6(C)) confirms the previous results. The ratio of participants whose responses are best described by a given model indicates that B-R describes the population at σ_g = 0.03 and at σ_π = 0.5 well. As in the bimodal dataset, MAP-R reflects the responses of some participants well at σ_g = 0.3. At σ_g = 0.1, the log likelihood performance (Fig 6A) of all models is similar but B-R_σ and MAP-R win by a small margin (see S1 Text for a subject-level analysis). Hence, B-R_σ and MAP-R perform best in the Bayesian random effects analysis (Fig 6C). No clear winner emerges at σ_g = 0.4.

Next, we revisit the question of whether participants sample or use a loss function. In the bimodal data, the loss function approach was at a disadvantage because it could only produce unimodal response distributions. This disadvantage is not present in the unimodal dataset. To make the performance of models in Fig 6(A) and 6(B) comparable, we use B-R with sampling as the baseline in both plots. Fig 6(B) shows that, averaging across participants and conditions, B-R with sampling outperforms the loss-based models. The large errors in (B) reflect large intersubject variability. The Bayesian random effects analysis in Fig 6(D) confirms that B-R also wins at the subject-level at σ_g = 0.03 and at σ_π = 0.5. One exception is L:B-R_σ at σ_g = 0.1. Indeed, the subject-level analysis (S1 Text) shows that in terms of the averaged log likelihood at middle and high noise σ_g ∈ (0.1, 0.4) participant 1 is an outlier. For other participants, the performance of B-R with sampling and the loss-based models is very similar. Despite the higher intersubject variability in the case of the unimodal dataset than in the bimodal dataset, the unimodal analysis provides convincing evidence of the superiority of B-R with sampling over other models considered here. In contrast to the model comparison in the bimodal analysis, in the unimodal case B-R clearly wins the model comparison at σ_π = 0.5.

Stimulus generation from the bimodal prior

Here, we describe in detail how stimuli are generated. On the j^th trial, participants are presented with a stimulus consisting of N = 4 points in a 2-dimensional space: . For the main experiment (with σ_π = 0.1, see below), we fixed the x-values to (−0.3, −0.1, 0.1, 0.3) respectively. For the additional experiment (with σ_π = 0.5, see below), we drew the x-values from Gaussians with means (−0.18, −0.09, 0, 0.09) and standard deviation 0.09 but resampled if the minimal distance was less than 0.1 between any two points. In both cases, we then generated the y-coordinates from a Gaussian generative model with a parabolic non-linearity and the generative parameter, w_j: (4) The parameter w_j is drawn from a mixed Gaussian prior (5) where the parameter set consists of the mean μ_π = 1, mixing coefficient c = 1/2 and the standard deviation σ_π = 0.1 for the main experiment and σ_π = 0.5 for the additional experiment. We denote the total set of hyperparameters (suppressed for notational clarity), from the prior and the generative probability, by . Each parameter w_j corresponds to a generative parabola. Given this model and given a stimulus D_j, we asked participants to predict the y-component y^⋆ at x^⋆ = 2, which is equivalent to mentally fitting a parabola to the four stimulus points and estimating the point of intersection with a vertical line at x^⋆.

To train participants on the generative model and the prior, we showed participants the generative parabola after each trial. In the main experiment, we showed a set of 20 unique stimuli for each of the three noise levels σ_g ∈ {0.03, 0.1, 0.4}, and each unique stimulus was repeated 20 times. We denote the set of the 20 responses to the j^th stimulus as . This amounts to a total of 400 trials per noise level. The order of the stimuli was randomized. For the additional experiment, we set σ_g = 0.1 and showed 10 unique stimuli 10 times. Fig 1 shows the experimental paradigm.

Regression models

In each trial, we model the participant’s computation by a consecutive inference and prediction step. During the inference step, the model assumes that the participant infers information about the quadratic parameter w_j based on the presented data (i.e., stimulus) D_j. The inferred information is then used for a subsequent prediction y^⋆. We describe the participant’s overall task as computing the predictive distribution: .

Prior regression (P-R) is our null model. P-R assumes that participants make predictions based on their prior belief but disregard information from the stimulus: (6) Maximum likelihood regression (ML-R) relies only on the likelihood maximizing parameter, w_ML: (7) Maximum a posteriori regression (MAP-R) uses the parameter that maximizes the posterior p(w|D_j) = p(D_j|w)π(w)/p(D_j): (8) Bayesian regression (B-R) uses the entire posterior for making predictions by marginalizing over it: (9) Bayesian regression with noise estimate (B-R_σ) loosens the assumption that participants treat σ_g as a hyperparameter and instead assumes they use an estimate on a trial-by-trial basis. Using the maximum likelihood estimator and the number of points M = 4: (10) After substituting the estimate for the hyperparameter σ_g in Eq (9), the posterior predictive distribution is computed analogously to B-R.

Participants’ internal noise

To predict the participants’ responses r from the regression models’ output y^⋆, we had to account for the internal noise of the participants. We did this by showing a noise-free stimulus 20 times and fitting a Gaussian with variance to each participant’s response distribution: . To be robust against outliers, we took the average of the 16% and 84% percentiles of the response distribution as motor noise. The predicted response distribution is then (11)

Model comparison

We used the log-likelihood and the variance to compare the predicted and empirical response distributions.

Determining the components of the variance in the response data

To compare the components of the predictive variance in Eq 3 to data, we make the following definitions for a set of response . The empirical mixing coefficient is the fraction of positive responses: If only one of the modes is present in the data (c ∈ {0, 1}) the bimodal contribution vanishes and we do not require the means for the total variance. If both modes are present we compute their means: We define the unimodal variance contribution as the variance of the dominant mode: where is the set of responses in the dominant mode, i.e., the mode that has the majority of responses. If no dominant mode exists we omit the stimulus. We did not use the inferior mode to have sufficient samples (at least 11) to estimate the variance.

The unimodal dataset

To obtain a unimodal dataset from the full dataset, we consider only responses and model predictions if they have the same dominant mode, i.e., parabolas facing either upwards or downwards. We define the dominant response mode as the one containing more than half of the responses and the dominant mode of the model as the one carrying more than half of the probability mass. For example, if 11 responses fall into the upper mode but the models predict a downward parabola, all responses are disregarded. However, if the models predict an upward parabola the 11 responses enter the unimodal dataset. Note that the symmetric prior ensures that B-R, B-R_σ and MAP-R share a dominant mode. Because the models use the same likelihood term, they process the stimulus as evidence for the same mode and break the symmetry in the same direction.

Averaged over participants, the fraction of trials per condition retained for the unimodal dataset is 0.994, 0.849, 0.596 and 0.703 for σ_g ∈ {0.03, 0.1, 0.4} with σ_π = 0.1 and σ_g = 0.1, σ_π = 0.5, respectively.

Participants

Seven naive participants (3 females, 4 males, ages 21-27) participated in the main experiment and four naive participants (all males, ages 21-30) took part in the second experiment. The experiments were programmed using custom software implemented in MATLAB. Stimuli were presented on a 1920x1080 (36 pixels/cm) monitor with a refresh rate of 120 Hz. Participants viewed the display binocularly. Each trial comprised a fixation dot presented for 1 s followed immediately by presentation of the stimulus (with 5 arcmin point diameter). Participants moved a red point up or down using the up and down arrow keys to indicate the vertical position of the parabola at the given horizontal location. See S1 Text for more details.

Ethics statement

All participants gave informed consent in accordance with protocol 384/2011 “Commission cantonale d’éthique de la recherche sur l’être humain”. Participants provided written consent prior to the experiment.