Bayesian inference in physics

Bayesian inference provides a consistent method for the extraction of information from physics experiments even in ill-conditioned circumstances. The approach provides a unified rationale for data analysis, which both justifies many of the commonly used analysis procedures and reveals some of the implicit underlying assumptions. This review summarizes the general ideas of the Bayesian probability theory with emphasis on the application to the evaluation of experimental data. As case studies for Bayesian parameter estimation techniques examples ranging from extra-solar planet detection to the deconvolution of the apparatus functions for improving the energy resolution and change point estimation in time series are discussed. Special attention is paid to the numerical techniques suited for Bayesian analysis, with a focus on recent developments of Markov chain Monte Carlo algorithms for high-dimensional integration problems. Bayesian model comparison, the quantitative ranking of models for the explanation of a given data set, is illustrated with examples collected from cosmology, mass spectroscopy, and surface physics, covering problems such as background subtraction and automated outlier detection. Additionally the Bayesian inference techniques for the design and optimization of future experiments are introduced. Experiments, instead of being merely passive recording devices, can now be designed to adapt to measured data and to change the measurement strategy on the fly to maximize the information of an experiment. The applied key concepts and necessary numerical tools which provide the means of designing such inference chains and the crucial aspects of data fusion are summarized and some of the expected implications are highlighted.

Figure 1

Three measurements of a physical quantity (e.g., mass) with Gaussian uncertainty, the allowed range being limited to [0; 0.2].

Figure 2

Posterior distribution $p(m|\mathbf{d},\sigma ,I)$ for three measurements with Gaussian uncertainty. The mean value $⟨m_0⟩$ is denoted by a vertical bar and the extent of the smallest 95% credibility region is indicated by the shaded area underneath the posterior distribution. Additionally the posterior distribution for the case of one measurement with a Cauchy (Lorentzian-) error distribution instead of a Gaussian distribution is given by a dashed line.

Figure 3

Left panel: The density of the slope is constant. This results in a nonuniform distribution for the angle between the straight lines and the $x$ axis. Right panel: The angular density is kept constant. Adapted from 115.

Figure 4

HD 73526 radial velocity measurement plotted from data given by 439.

Figure 5

HD 73526 radial velocity measurement data superimposed with the best-fit model radial velocity. Adapted from 191.

Figure 6

Observations of the onset of cherry blossoms at Geisenheim, Germany from 1900 to 2002 in terms of days after the beginning of the year. Two observations (1946 and 1949) are missing. Adapted from 119.

Figure 7

Data of cherry blossom onset with normalized change point probability (full squares) superimposed (see right vertical scale). Adapted from 119.

Figure 8

Posterior density for the signal rate $s$ for $N_{\mathrm{on}}=3$, $N_{\mathrm{off}}=7$, and $t=1$ s and uniform priors. The maximum is for $s=0$, corresponding to a vanishing signal rate, but the gray shaded 95% credible interval extends up to $s=4{\mathrm{s}}^{-1}$.

Figure 9

RBS spectra of thin Co and Cu films on a Si substrate, measured with 2.6 MeV He. The Cu spectrum (right-hand side) is deconvolved with the apparatus function obtained from the Co spectrum (left-hand side). The two Cu isotopes are clearly resolved with measured abundances close to the natural abundances. Adapted from 147.

Figure 10

RBS data of the sample before and after plasma exposition. The signal edge position at around 550 keV is shifted toward lower energies due to the plasma exposition. At the same time the peak intensity increases relative to the bulk signal. From 455.

Figure 11

Reconstructed depth profiles and asymmetric confidence intervals from the RBS spectra shown in Fig. 10. (a) The sample composition before exposure (lines are added to guide the eye), and (b) the sample composition after exposure. From 459.

Figure 12

Example of a two-dimensional probability distribution with properties complicating the sampling: More than one maximum, regions of high probability are separated, non-Gaussian shape, no alignment with the coordinate system.

Figure 13

Rejection sampling: First, sample a candidate $x^{(0)}$ from $p(x)$ and a uniform variable $u$. Then accept this candidate if $t(x^{(0)})\ge ukp(x^{(0)})$ (as in the situation shown here); otherwise reject the candidate.

Figure 14

Sample of 2000 independent samples drawn by rejection sampling from the 2D test distribution shown in Fig. 12.

Figure 15

Example of a typical MCMC sample, showing indication of random-walk behavior. Solid lines indicate accepted proposals, dashed lines rejected proposals. The underlying potential is the one displayed in Fig. 12, upper right part, and represented by ellipsoidal contour lines.

Figure 16

Influence of the width $\sigma$ of a Gaussian proposal function on the autocorrelation of the samples. Too small and too large widths result in a very slow decay of the autocorrelation. If the width of the proposal function is too small then the chain hardly explores the parameter space. If instead the width of the proposal function is too large then almost all proposed changes are rejected again causing insufficient exploration.

Figure 17

Influence of the width of the proposal function on the burn-in and convergence of the Markov chains. The simulation with a narrow proposal function (solid black line) exhibits a larger autocorrelation length of an arbitrarily chosen parameter $x_1$ and a slower convergence to the stationary distribution compared to the simulation with a 5 times larger proposal width.

Figure 18

Gibbs sampling of a correlated two-dimensional probability distribution. The alternating updates of the two variables are responsible for the “rectangular” movement of the chain. The step size is governed by the respective conditional distributions. Unlike in standard-MCMC methods all proposals are accepted. The underlying potential is the one displayed in Fig. 12, upper right part, and represented by ellipsoidal contour lines.

Figure 19

Example of a typical parallel tempering run with five replicas at different temperatures $0={\beta }_0<{\beta }_1<\cdots <{\beta }_4=1$. The average time of the replicas to swap from $\beta =0$ to $\beta =1$ and back is a measure of the efficiency of a parallel tempering scheme.

Figure 20

Schematics of one-dimensional slice sampling: For a given location $x^{(t)}$ a uniform random number $u$ is sampled from [0, $p(x^{(t)})$] and an interval [$x_{\min }$, $x_{\max }$] is determined such that at both locations $\{x_{\min },x_{\max }\}$ the probability distribution is smaller than $u$. In this interval a position $x^{\prime }$ is accepted as a new location if $p(x^{\prime })>u$.

Figure 21

Example of a two-dimensional Hilbert curve on a ${32}^2$ grid. The space-filling curve can be used to map distributions into a lower-dimensional space with arbitrary precision.

Figure 22

Part of the two-dimensional probability distribution of Fig. 12 being mapped by a ${256}^2$ Hilbert curve (similar to the one in Fig. 21) into one dimension. The index $i$ labels the points of the discretized 2D Hilbert curve in consecutive order, and $p(i)$ gives the probability at the corresponding location [$x_1(i)$, $x_2(i)$]. The mapping results in a ragged shape of the distribution; however, this is no obstacle for most MCMC algorithms.

Figure 23

Illustration of the dimension matching in RJMCMC. The lower-dimensional parameter spaces are padded with additional dimensions to allow a comparison with the higher-dimensional models.

Figure 24

Example of the two different cases for expectation value computation as discussed. (a) The probability distribution $p(x)$ displays more structure than the function $f(x)$ of which the expectation value is taken. This is the easier case. (b) The expectation value is computed from a function which is not well matched to the probability distribution, often the case in the evaluation of the marginal likelihood, requiring the use of advanced methods (see Sec. 4e3).

Figure 25

Nested sampling: (a) The sample points are uniformly sampled from the prior distribution. (b) The sample points are sorted with respect to prior volume enclosing likelihood $\mathrm{mass}\le L(\xi )$.

Figure 26

Sky map of the cosmic microwave temperature fluctuations around the mean value of 2.725 K in galactic coordinates (Mollweide projection) after subtraction of foreground sources. The foreground-reduced internal linear combination map is based on the 5-yr WMAP data. From 209.

Figure 27

Linear interpolated reconstruction of the primordial curvature fluctuation spectrum. The bins are logarithmically spaced. The errors show the 68% and 95% constraints and the black diamonds show the mode of the likelihood. The dashed vertical line on the left-hand side shows the values for $k=0$. Only the amplitude was allowed to vary at each of the nodes. From 413.

Figure 28

Low-$l$ multipoles and $1\sigma$ error bars from three releases of WMAP$WMAP$ data. The best-fitting fiducial power spectrum based on WMAP5 inferences is indicated together with the associated cosmic variance limits. The multipol values are slightly shifted for clarity. From 53.

Figure 29

Piecewise linear interpolated reconstruction of the primordial spectrum with different degrees of flexibility. Only the amplitude was allowed to vary at each of the support points (indicated by filled symbols). In model 1 only the amplitude of a flat spectrum could be varied (shown as a dashed line in both graphs), in model 2 the slope could also be adjusted (two support points, solid line, upper graph), and model 3 corresponds to a piecewise linear model with two segments (corresponding to three support points, solid line, lower graph). The Bayes factors are $B_{21}=2$ and $B_{31}=3$ in favor of the more complex models compared to the flat model 1. Adapted from 53.

Figure 30

Comparison of data computed by two different models (but with the same number of radicals) with measured mass spectrometer data. Model 6 additionally incorporates the species ${\mathrm{C}}_4{\mathrm{H}}_2$ and ${\mathrm{C}}_4{\mathrm{H}}_6$ compared to model 1 which contains only up to ${\mathrm{C}}_3{\mathrm{H}}_x$ species. The former model provides a near perfect fit to the measured data even at high masses.

Figure 31

The natural logarithm of $p(E|\mathbf{D},\mathbf{S},I)$ is displayed using full dots; the scale is given on the left ordinate (indicated by the upper arrow). The misfit between data and model for different combinations of six free radicals (${\mathrm{C}}_2{\mathrm{H}}_5$, ${\mathrm{CH}}_3$, H, ${\mathrm{C}}_2{\mathrm{H}}_3$, CH, and ${\mathrm{CH}}_2$) (while keeping the set of nonradical species constant) is indicated with open symbols, the corresponding scale is given on the right ordinate as the averaged ${\chi }^2$ normalized by the number of data points, indicated by the lower arrow. The abscissa shows the number of radicals involved in the model, which were taken in the given order. Lines are to guide the eye. From 244.

Figure 32

Flux dependence of the chemical erosion yield of graphite under hydrogen irradiation. The abscissa shows the flux of deuterium, and the ordinate shows the carbon erosion yield in eroded carbon atoms per impinging deuterium ion. The data set $\delta$ is represented by circles. Filled circles correspond to the subset for which the fitting curve (solid line) is valid ($E_0=30\mathrm{eV}$, $T=600\mathrm{K}$). Open triangles and squares represent data from the experiments PSI-I (195) and PISCES (464), respectively, while the full symbols show the data sets after multiplication with the corresponding scale factors (0.72 and 0.32). Error bars show the assigned experimental error. The gray shaded area is the confidence range. From (360.

Figure 33

(a) An $MVV$ Auger spectrum for iron. The estimated background shown is obtained for the transformed spectrum shown in (b). A logarithmic transformation of the Auger spectrum reduces the curvature of the background. The estimated background is shown as a solid line. The eight support points of the spline are indicated by filled circles. (c) The signal obtained by subtracting the data and the background. A secondary peak is present at an energy of 86, 39 eV above the $M_{2,3}VV$ Auger transition, substantiated by the autocorrelation of the signal vs energy difference (see inset). Adapted from 145.

Figure 34

Traditional approach to data fusion; based on some initial guess of the required parameters the respective data of each diagnostic are evaluated to yield individual best-fit parameters. These are combined, and additional constraints (e.g., positive density) are taken into account and used as starting parameters for the next iteration cycle. This process is repeated until a self-consistent solution is obtained.

Figure 35

Bayesian approach to data fusion; based on the prior distribution of the required parameters the respective likelihood distribution of each diagnostic is evaluated. Then, using Bayes theorem the posterior distribution of the parameters is computed, automatically taking into account all interdependencies of the diagnostics.

Figure 36

Posterior distribution of electron density at the position $z=6\mathrm{cm}$ for a plasma discharge derived from combined evaluation of Thomson scattering data and x-ray diagnostic.

Figure 37

Posterior distribution of electron temperature and electron density at the position $z=6\mathrm{cm}$ for a plasma discharge derived from Thomson scattering data.

Figure 38

Posterior distribution of electron temperature given by a soft-x-ray diagnostic. The constant values along the vertical direction indicate that the measurement provides no information about the electron density $n_e$.

Figure 39

Two-dimensional posterior distribution of electron temperature and electron density at the position $z=6\mathrm{cm}$ for a plasma discharge derived from Thomson scattering data.

Figure 40

Representation of a Bayesian network as a directed acyclic graph (DAG), the causal directions indicated by arrows.

Figure 41

Uncertainty of the magnetic flux surfaces represented by posterior samples. Adapted from 111.

Figure 42

Internal representation of the trajectory (continuous line) derived from odometer data of a robot exploring a building superimposed onto a map of the building (dark areas represent walls). The starting point is in the upper right corridor of the building, indicated by an open circle. The mismatch is continuously increasing due to accumulation of orientation angle uncertainties. Adapted from 151.

Figure 43

2D projection of the 3D-probability distribution visualizing the accumulated localization uncertainty after moving from a precisely known initial position with several 90° turns. The robot is more likely to be in darker areas. The line describes the ideal trajectory of the robot. Adapted from 151.

Figure 44

The extended Kalman filter (EKF) applied to the robot mapping problem. The robot’s path is indicated by the dotted line, and its estimations of its own position are shown as shaded ellipses which may deviate from the true position (indicated by crosses). Eight distinguishable landmarks of unknown location are indicated by small squares, and their location estimates by the robot are shown as open ellipses. During the path the robot’s positional uncertainty is increasing, as is its uncertainty about the landmarks it encounters. Adapted from 387.

Figure 45

Once the robot senses again the first landmark (indicated by an arrow) the uncertainty of all landmark estimates decreases (indicated by the reduced size of the unfilled ellipses), thanks to the fact that all the estimates are correlated. In addition the new position estimate of the robot itself is also much more precise. Adapted from 387.

Figure 46

The true velocity curve used for the generation of mock data is displayed as a dashed line. From this velocity curve ten simulated observations distorted with Gaussian noise were generated. These observations are given as points with error bars. From 286.

Figure 47

Posterior samples from the probability distribution for two velocity curve parameters, the spread indicating the uncertainty in the marginalized parameters. The time for an orbital period is given by $\tau$. The orbital eccentricity $e$ describes the amount by which an orbit deviates from a perfect circle: $e=0$ is a perfectly circular orbit and $e=1$ corresponds to an open parabolic orbit. From 286.

Figure 48

The true velocity curve used for the simulation of the data is given as a dashed line. A posterior sample of 15 predicted velocity curves based on the observations is given by thin solid lines; their spread is a measure of the uncertainty of the predicted velocity. The expected information gain for a further measurement at each time is indicated by a thick solid curve; the information gain is given on the right axis.Note that the positions of largest uncertainty and largest information gain coincide. From 286.

Figure 49

Samples from the probability distribution for two velocity curve parameters after one additional measurement at the best time. Comparison with Fig. 47 reveals that the parameters have a significantly increased precision, the samples are less disperse. The time for an orbital period is given by $\tau$. The orbital eccentricity $e$ describes the amount by which an orbit deviates from a perfect circle: $e=0$ is a perfectly circular orbit and $e=1$ corresponds to an open parabolic orbit. From 286.

Figure 50

Simulated yield data of a $\mathrm{D}({}^3\mathrm{He},p){}^4\mathrm{He}$ nuclear reaction analysis of a tungsten sample with an exponentially decaying deuterium concentration profile. The varying intensity reflects the interplay of an energy dependent cross section, ion beam range, and concentration profile. Adapted from 460.

Figure 51

Three cycles of the experimental design process. On the left-hand side 1000 samples drawn from the posterior distribution $p(a_1,a_2|\mathbf{d},\eta )$ are displayed. On the right-hand side the expected utility is plotted and the maximum is indicated by a circle. The corresponding abscissa value is the suggested next measurement energy. Performing that measurement yields the posterior samples (black dots) given in the second row, left-hand side. The previous posterior samples are given in gray for comparison. On the right-hand side the new expected utility as a function of beam energy is given. Adding a new measurement with the proposed energy results in the posterior sample displayed in the lower left panel. The previous posterior samples are also given in different colors. The posterior volume has already decreased significantly. Therefore the best EU for the next measurement is lower than before. Adapted from 460.

Figure 52

(A), (C) The area to be searched together with the position of the white disk to be detected. In (A) a first measurement has been taken (small filled black circle in the upper right). Since at that location background only was detected, a number of possible positions of the white disk can immediately be excluded. This is visualized by 150 circles sampled from the updated posterior. These represent possible disk positions and sizes consistent with the measurement(s). (C) The situation after several measurements: The algorithm is getting close to a solution. The possible circle parameters are already narrowed down. The scatter of possible disk positions is already pretty small. (B) and (D) The selection algorithm for the location of the next measurement. The panels display the information gain for each measurement location (dark colors indicate uninformative locations) and the location with the highest information gain is selected for the next measurement (indicated by arrows). Further details are given in the text. From 256.

Figure 53

Parameters of interest for the design of a multichannel interferometer for W7-X. The parameters are varied according to different high confinement regimes. The different lines (solid, dashed, and dotted) represent different realizations of physical scenarios On the left-hand side the maximum density of the plasma is varied, keeping the edge position constant. On the right-hand side the maximum density and the maximum slope are kept constant, varying the position of the plasma edge only. Adapted from 127.

Figure 54

Design result for four interferometry chords with respect to the measurement of high confinement regimes: Four beam lines in the interferometry plane (left) and optimal four beam configuration without technical constraints (right) superimposed onto a cross section of the W7-X plasma. Adapted from 128.