Model and Statistical Framework

This page is for readers who want to understand the statistical machinery behind bayesian_listener. Equation numbers refer to Barumerli et al.[1] and Barumerli et al.[2]; the two papers use consistent notation.

Model pipeline

The model simulates a single static-sound-localisation trial in three stages:

1. Feature extraction. The binaural stimulus generated by a source at direction \(\boldsymbol{\varphi}\) is characterised by a vector of noisy spatial features (Eq. 1 of Barumerli et al.[1]):

   \[\mathbf{t} = [x_\mathrm{itd},\, x_\mathrm{ild},\, \mathbf{x}_{L,\mathrm{mon}},\, \mathbf{x}_{R,\mathrm{mon}}] + \boldsymbol{\delta}\]

   where \(x_\mathrm{itd}\) is the interaural time difference (ITD), \(x_\mathrm{ild}\) the interaural level difference (ILD), and \(\mathbf{x}_{L/R,\mathrm{mon}}\) the monaural log-amplitude spectra for the left and right ears. The additive noise \(\boldsymbol{\delta} \sim \mathcal{N}(\mathbf{0}, \boldsymbol{\Sigma})\) represents perceptual uncertainty, with diagonal covariance (Eq. 2 of Barumerli et al.[1]):

   \[\begin{split}\boldsymbol{\Sigma} = \begin{pmatrix} \sigma_\mathrm{itd}^2 & 0 & 0 \\ 0 & \sigma_\mathrm{ild}^2 & 0 \\ 0 & 0 & \sigma_\mathrm{mon}^2 \mathbf{I} \end{pmatrix}\end{split}\]

2. Bayesian inference. The listener’s internal estimate of the source direction is the maximum a-posteriori (MAP) direction \(\hat{\varphi}'\) obtained by combining the sensory likelihood with an elevation prior (Eqs. 4–6 of Barumerli et al.[1]):

   \[p(\mathbf{t} \mid \boldsymbol{\varphi}) = \mathcal{N}(\mathbf{t} \mid \mathbf{s}(\boldsymbol{\varphi}),\, \boldsymbol{\Sigma})\]
   \[p(\boldsymbol{\varphi}) \propto \exp\!\left( -\frac{\epsilon^2}{2\,\sigma_\mathrm{prior}^2}\right)\]
   \[\hat{\varphi}' = \arg\max_{\boldsymbol{\varphi}}\; p(\mathbf{t} \mid \mathbf{s}(\boldsymbol{\varphi}))\, p(\boldsymbol{\varphi})\]

   where \(\epsilon\) is the elevation angle of \(\boldsymbol{\varphi}\) and \(\mathbf{s}(\boldsymbol{\varphi})\) is the noiseless template at that direction, interpolated from the listener’s HRTF.

3. Motor noise. The final pointing response \(\hat{\varphi}\) is the MAP estimate perturbed by direction-independent motor noise (Eq. 7 of Barumerli et al.[1]):

   \[\hat{\varphi} = \hat{\varphi}' + \mathbf{m}, \quad \mathbf{m} \sim \mathrm{vMF}(\mathbf{0},\, \kappa_m)\]

   where \(\kappa_m\) is the von Mises–Fisher concentration parameter, related to the motor-noise standard deviation \(\sigma_m\) in degrees via \(R = I_1(\kappa_m)/I_0(\kappa_m) = \exp(-\sigma_m^2/2)\).

Likelihood function

For a localisation experiment with \(R\) trials, the model likelihood is (Eq. 8 of Barumerli et al.[2]):

\[\mathcal{L}(\boldsymbol{\theta} \mid [\hat{\varphi}_1^*, \ldots, \hat{\varphi}_R^*],\, [\boldsymbol{\varphi}_1, \ldots, \boldsymbol{\varphi}_R]) = \prod_{i=1}^{R}\, p(\hat{\varphi}_i^* \mid \boldsymbol{\varphi}_i,\, \boldsymbol{\theta})\]

where \(\boldsymbol{\theta} = \{\sigma_\mathrm{itd}, \sigma_\mathrm{ild}, \sigma_\mathrm{mon}, \sigma_\mathrm{prior}, \kappa_m\}\).

Because the integral over internal estimates has no closed form (Eq. 9 of Barumerli et al.[2]), it is approximated via Monte Carlo with \(M\) samples (Eq. 10 of Barumerli et al.[2]):

\[p(\hat{\varphi}^* \mid \boldsymbol{\varphi},\, \boldsymbol{\theta}) \approx \frac{1}{M} \sum_{m=1}^{M} \mathrm{vMF}(\hat{\varphi}^* \mid \hat{\varphi}_m,\, \kappa_m)\]

\(M = 200\) samples is sufficient for stable likelihood approximation (see Fig. S2 of Barumerli et al.[2]).

Model comparison uses the Bayesian Information Criterion (Eq. 14 of Barumerli et al.[2]):

\[\mathrm{BIC} = k \ln(R) - 2\ln\mathcal{L}(\hat{\boldsymbol{\theta}})\]

where \(k\) is the number of free parameters.

Datasets

The package has been validated on the SONICOM Multi-Experiment Auditory Localisation Dataset. It contains localisation trials from 34 participants tested across seven HRTF conditions (individual measured, synthetic individual, best-match, worst-match, KEMAR, and random non-individual HRTFs) collected at the Turret Lab, Imperial College London. Responses are provided as CSV with per-trial azimuth, elevation, and great-circle error in both spherical and horizontal-polar coordinates. The validation and fitting of the model has been done on the individual measured condition.

Noise parameters

The model has five noise parameters. Two are fixed at literature values; three are estimated by the two-stage fitting procedure.

Symbol	Physical interpretation	Typical range	Identifiability	Fixed or fitted
\(\sigma_\mathrm{itd}\)	ITD perceptual noise (dimensionless)	— (fixed)	—	Fixed at 0.569
\(\sigma_\mathrm{ild}\)	ILD perceptual noise (dB)	— (fixed)	—	Fixed at 1.0 dB
\(\sigma_m\)	Motor noise standard deviation (degrees)	2°–36° (median ≈ 10°)	Good (\(r = 0.97\))	Fitted (stage 1, lateral only)
\(\sigma_\mathrm{mon}\)	Monaural spectral noise (dB)	2–15 dB	Moderate (\(r = 0.85\))	Fitted (stage 2, full sphere)
\(\sigma_\mathrm{prior}\)	Elevation prior width (degrees)	5°–90°	Poor–Moderate (\(r = 0.84\), positive bias)	Fitted (stage 2, full sphere)

Identifiability ratings and ranges are from the parameter recovery analysis in Barumerli et al.[2].

Template interpolation

The template \(\mathbf{s}(\boldsymbol{\varphi})\) is computed by interpolating listener-measured HRTF features onto a quasi-uniform spherical grid of \(T = 2{,}112\) directions (4° average spacing). Four methods are available:

• SHMAX (recommended): regularised spherical-harmonic (SH) interpolation with order selected per dataset to maximise numerical stability.

• SH: SH interpolation with a fixed order; can undersmooth high-frequency spectral detail.

• barycentric: barycentric interpolation on the Delaunay triangulation of the sampling grid; computationally lighter than SH.

• barumerli2023: the original MATLAB implementation; retains features only above the minimum measurement elevation and is included for backward compatibility.

Full-sphere spatial coverage and high-frequency spectral fidelity are the primary determinants of template quality; among full-sphere methods the choice of algorithm is secondary for dense measurement grids such as SONICOM (Barumerli et al.[2], Sec. 3.2).

Known limitations

The following limitations affect the interpretation of fitted parameters and model predictions.

Lateral accuracy bias.

The model cannot reproduce a non-zero mean lateral error by construction (the response noise is zero-mean). Observed lateral biases likely reflect pointing-apparatus calibration or individual motor asymmetries not captured by the current formulation.

Poor identifiability of \(\sigma_\mathrm{prior}\) at large values.

When \(\sigma_\mathrm{prior}\) exceeds roughly 70°, the prior becomes approximately uniform over elevation and its contribution to the posterior is indistinguishable from a flat prior. Recovered values show a systematic positive bias at large ground-truth values (mean bias +19.9°, Barumerli et al.[2]).

Motor noise degrades for near-chance polar responders.

\(\sigma_m\) is estimated from lateral responses; participants with near-chance polar performance carry little information about \(\sigma_\mathrm{mon}\), resulting in noisier spectral-noise estimates.

Variable trial counts across participants.

The behavioural dataset pools three experiments with different repetitions per direction (3, 6, or 9). Raw likelihood comparisons across participants should account for trial count via BIC.

Extrapolation below the measurement boundary (barumerli2023 only).

The barumerli2023 method assigns zero template mass below the lowest measured HRTF elevation, biasing the posterior for sources near or below the horizontal plane. Full-sphere methods avoid this by extrapolating features into the lower hemisphere.

Further reading

The standard localisation metrics implemented in bayesian_listener.metrics follow the definitions of Middlebrooks[3]. The gammatone filterbank used in compute_features follows the ERB-rate scale of Glasberg and Moore[4]. The barycentric (VBAP) interpolation in vbap_interpolate follows Pulkki[5], and parameter optimisation throughout bayesian_listener.fitting is performed with BADS Acerbi and Ma[6].

References

[1] (1,2,3,4,5)

R. Barumerli, P. Majdak, M. Geronazzo, D. Meijer, F. Avanzini, and R. Baumgartner. A Bayesian model for human directional localization of broadband static sound sources. Acta Acustica, 7:12, 2023. doi:10.1051/aacus/2023006.

[2] (1,2,3,4,5,6,7,8,9)

R. Barumerli, F. Brinkmann, E. Zanoni, A. Hoyer, L. Picinali, and M. Geronazzo. Statistical validation and full-sphere extension of a Bayesian model for human static sound localisation. Submitted to Acta Acustica, 2026. URL: https://arxiv.org/abs/2606.24367.

[3]

J. C. Middlebrooks. Virtual localization improved by scaling nonindividualized external-ear transfer functions in frequency. J. Acoust. Soc. Am., 106(3):1493–1510, 1999. doi:10.1121/1.427147.

[4]

B. R. Glasberg and B. C. J. Moore. Derivation of auditory filter shapes from notched-noise data. Hearing Research, 47(1–2):103–138, 1990. doi:10.1016/0378-5955(90)90170-T.

[5]

V. Pulkki. Virtual sound source positioning using vector base amplitude panning. J. Audio Eng. Soc., 45(6):456–466, 1997.

[6]

L. Acerbi and W. J. Ma. Practical Bayesian optimization for model fitting with Bayesian adaptive direct search. In Advances in Neural Information Processing Systems, volume 30. 2017.