Solar Bayesian Analysis Toolkit—A New Markov Chain Monte Carlo IDL Code for Bayesian Parameter Inference

Once the most credible value θ_MAP of the model parameters is determined, one can calculate the predictive distribution of observational data points (i.e., what the next observation D_new could be):

$$\begin{eqnarray}&&P({D}_{\mathrm{new}}| {\theta }_{\mathrm{MAP}})=P({D}_{\mathrm{new}}| \theta ={\theta }_{\mathrm{MAP}}).\end{eqnarray} \tag{ 4 }$$

However, Equation (4) does not account for the estimate θ_MAP being uncertain itself. This uncertainty comes from the observational errors and model limitations, and is the width of the posterior PDF in the vicinity of its global maximum. To account for all uncertainties correctly, the posterior predictive distribution

$$\begin{eqnarray}&&P({D}_{\mathrm{new}}| D)=\int P({D}_{\mathrm{new}}| \theta )P(\theta | D)d\theta \end{eqnarray} \tag{ 5 }$$

is used. It is usually broader than the distribution given by Equation (4) because of the additional uncertainties in θ.

The posterior predictive distribution can be used for two purposes. The first is to forecast future observations and to provide reliable prediction intervals, if the model allows for extrapolation in time. The second application is a so-called posterior predictive check, which allows for assessment of the consistency of the chosen model with the observations in terms of confinement of the model in the data space. A reliable model should produce a narrow distribution predicting possible observations of the same process to be close to the actual data points.

Before discussing model comparison criteria let us list properties of a good model:t

• 1.  
  the data points are well approximated by the best fit provided by a model;
• 2.  
  the model is sufficiently simple and does not "overfit" the data;
• 3.  
  the model is well confined in the parameter space and provides sufficiently small uncertainties for free parameter estimates;
• 4.  
  the model is well confined in the data space and does not predict possible observations far away from the actual data points.

The meaning of the words "well" and "sufficient" is, of course, subjective. However, if one has to select the best model from several available options, these criteria become objective and can be quantitatively assessed.

A problem of model comparison often appears in a data analysis workflow. This problem is traditionally used by applying one of the best-fit comparison criteria such as the classical (reduced) χ² and likelihood ratio tests or the more elaborated Akaike information criterion (AIC) and Bayesian information criterion (BIC). The latter two are based on information theory and Bayesian statistics, respectively. All criteria mentioned above are based on comparison of the deviation of the found fit from the data points. The AIC, BIC, and reduced χ² have a special term to penalize models with larger number of free parameters.

Despite the great success, all best-fit comparison criteria, including the BIC and AIC, have a major drawback. They compare only the best fits, while the model complexity is accounted for only partially by having an additional term dependent on the number of free parameters. Thus, the best-fit comparison criteria address only point 1 and partially point 2 in the list above.

The Bayesian approach provides us with a quantitative model comparison technique which transparently addresses all properties of a "good" model mentioned above. The technique is based on calculating the Bayesian evidence (2).

To understand the meaning of Bayesian evidence, let us use the Bayes theorem to compute the posterior probability of a model M_i in the case of N competing models M₁, M₂ ... M_N:

$$\begin{eqnarray}&&P({M}_{i}| D)=\displaystyle \frac{P(D| {M}_{i})P({M}_{i})}{P(D)},\end{eqnarray} \tag{ 6 }$$

where P(M_i) is the prior probability of model M_i. The likelihood term P(D∣M_i) is nothing but Bayesian evidence for model M_i defined as in (2). P(D) is a normalization constant and can be computed as

$$\begin{eqnarray}&&P(D)=\displaystyle \sum _{i=1}^{N}P(D| {M}_{i})P({M}_{i}).\end{eqnarray} \tag{ 7 }$$

After applying the Bayes theorem, the Bayesian evidence becomes a posterior probability for a model to be true under condition of observed data D.

Thus, posterior model probabilities (6) are a complete solution for the model comparison within the Bayesian approach. It transparently addresses all points listed in the beginning of this section, accounts for prior model probabilities and determines the posterior probabilities of every competing model. These probabilities allow for the quantitative model comparison and should be interpreted in a Bayesian sense as a degree of belief, where 0 is impossible and 1 is absolutely true.

The key ingredient in the Bayesian model comparison is the Bayesian evidence. Unlike best-fit metrics (χ², likelihood maximum value, AIC, BIC, and others), it is an integral over the whole parameter space and, thus, is rather a function of the whole model than a best-fit metric.

Quantitative comparison of two models M₁ and M₂ with equal prior probabilities can be done by calculating the Bayes factor (Jeffreys 1961) which is the ratio of Bayesian evidences (and posterior probabilities) for models M₁ and M₂:

$$\begin{eqnarray}&&{B}_{12}=\displaystyle \frac{P(D| {M}_{1})}{P(D| {M}_{2})},\end{eqnarray} \tag{ 8 }$$

where the evidences P(D∣M₁) and P(D∣M₁) are calculated according to Equation (2). Traditionally, the doubled natural logarithm of this factor is used, i.e.

$$\begin{eqnarray}&&{K}_{12}=2\mathrm{ln}{B}_{12},\end{eqnarray} \tag{ 9 }$$

where values of K₁₂ greater than 2, 6, and 10 correspond to "positive", "strong," and "very strong" evidence for model M₁ over model M₂, respectively (Kass & Raftery 1995). The calculation of corresponding posterior model probabilities using (6) and assuming P(M₁) = P(M₂) = 0.5 gives P(M₁∣D) of 73%, 95%, and 99% for "positive," "strong," and "very strong" evidence, respectively.

To generate a large number of samples from the posterior distribution, SoBAT uses the MCMC technique. The marginalized posterior PDFs are then approximated by the histograms of these samples.

The MCMC sampling algorithm is the most important part of our code. It can generate samples from the posterior distribution using any target function f(θ) which is proportional to the posterior PDF P(θ∣D) and is a known continuous function that can be calculated for any value of θ. Thus, the knowledge of the normalization constant (Equation (2)) is not required for the inference.

Our sampling algorithm is the classical random walk Metropolis–Hasting sampler with the multivariate normal distribution used as a proposal distribution. Its covariance matrix $\hat{\sigma }$ is automatically tuned to keep the acceptance rate in the range of 10%–90% during the whole sampling procedure. In order to generate the whole sequence of samples (chain) with the same proposal distribution, we restart the sampling procedure every time when the proposal distribution is tuned. The detailed description of the algorithm is given below.

• 1.  
  Initialize the starting point in the parameter space, Θ₀.
• 2.  
  Estimate the local covariance matrix $\hat{\sigma }$ for θ = Θ₀.
• 3.  
  Simulate the proposed sample Ξ_i from the multivariate normal distribution $N({{\rm{\Theta }}}_{i},\hat{\sigma })$ with the expected value Θ_i and covariance matrix $\hat{\sigma }$.
• 4.  
  Compute the ratio R = f(Ξ_i∣D)/f(Θ_i∣D).
• 5.  
  Pick a random number ε between 0 and 1.
• 6.  
  Produce a new sample Θ_i+1:

  $$\begin{eqnarray*}\left\{\begin{array}{ll}\mathrm{accept}:\ {{\rm{\Theta }}}_{i+1}={{\rm{\Xi }}}_{i};{N}_{{\rm{a}}}={N}_{{\rm{a}}}+1; & (\mathrm{if}\ \varepsilon \leqslant R)\\ \mathrm{reject}:\ {{\rm{\Theta }}}_{i+1}={{\rm{\Theta }}}_{i};{N}_{{\rm{r}}}={N}_{{\rm{r}}}+1; & (\mathrm{if}\ \varepsilon \gt R).\end{array}\right.\end{eqnarray*}$$
• 7.  
  Calculate the acceptance rate r = N_a/(N_r + N_a).
• 8.  
  If r < 10% or r > 90%⁵ then set Θ₀ = Θ_i+1 and go to step 2.
• 9.  
  Repeat steps 3–8 until the desired number of samples is generated.
• 10.  
  Return all collected samples Θ_i as a result.

After several restarts, the sampling algorithm usually finds the maximum probability area and stabilizes there with the acceptance rate at about 10%–90%. Note that there is no guarantee that the algorithm will find the global maximum for a given number of iterations. Therefore, we recommend providing a rather good initial guess and generating a sufficiently large number of samples.

3.1.1. Burning in Stage

The selection of the proposal distribution is essential for constructing an effective Metropolis–Hastings sampler. The developed code uses the multivariate normal distribution with the expected value μ = Θ₀ and the covariance matrix $\hat{\sigma }$, which is tuned to reflect the local properties of the parameter space and to achieve an optimal acceptance rate. The algorithm of the calculation of the optimal covariance matrix $\hat{\sigma }$ is given below.

• 1.  
  Initialize variables.

  • (a)  
    Θ₀—a position in the parameter space
  • (b)  
    $\hat{\sigma }$—an initial guess for the covariance matrix
  • (c)  
    S—an array to store generated samples
• 2.  
  Simulate the proposed sample Ξ_i from the multivariate normal distribution $N({{\rm{\Theta }}}_{0},\hat{\sigma })$ with the expected value Θ₀ and covariance matrix $\hat{\sigma }$.
• 3.  
  Compute the ratio $R=\min \left(\tfrac{f({{\rm{\Xi }}}_{i}| D)}{f({{\rm{\Theta }}}_{0}| D)},\tfrac{f({{\rm{\Theta }}}_{0}| D)}{f({{\rm{\Xi }}}_{i}| D)}\right)$.
• 4.  
  Generate a random number ε between 0 and 1.
• 5.  
  If ε ≤ R, accept and save sample S ← Ξ_i; N_a = N_a + 1 or reject it N_r = N_r + 1 otherwise.
• 6.  
  Calculate the acceptance rate r = N_a/(N_r + N_a).
• 7.  
  Tune $\hat{\sigma }$ for better acceptance rate

  • (a)  
    if r = 0 during 500 subsequent iterations, set $\hat{\sigma }=0.5\hat{\sigma }$
  • (b)  
    if r > 50%, set $\hat{\sigma }=1.1\hat{\sigma }$
• 8.  
  If more than 500 samples were accepted, set $\hat{\sigma }\,=\mathrm{covariance}(S)$
• 9.  
  Repeat steps 2–8 until the desired number of samples is generated.
• 10.  
  Return covariance(S) as a result.

The code allows evidences to be calculated by numerical evaluation of the integral given by Equation (2). The ratio of evidences for two models is the Bayes factor and can be interpreted as described in Section 2.2. The numerical integration of Equation (2) is implemented using the importance sampling Monte-Carlo technique (Hastings 1970). As an importance function, we use a multivariate Gaussian with the covariance matrix computed from the simulated MCMC samples from the posterior distribution.

To compute evidence for a given model, SoBAT offers the MCMC_EVIDENCE function. The function has three required parameters:

• 1.  
  f(θ)—a function computing the natural logarithm of a target function proportional to the posterior PDF;
• 2.  
  S_i, i = 1..N_s—samples simulated from the posterior by the MCMC_SAMPLE function;
• 3.  
  N—number of iterations for the Monte-Carlo integration.

The importance sampling Monte-Carlo integration is interpreted in the following form:

• 1.  
  Estimate the covariance matrix $[\hat{\sigma }]$ and the expected value [μ] from the posterior samples. The PDF n(θ) of the the multivariate normal distribution ${ \mathcal N }\left(\mu ,\hat{\sigma }\right)$ will be used as the importance function.
• 2.  
  Repeat N times (i = 1..N):

  • (a)  
    simulate a position⁶ θ_i in the parameter space from the multivariate normal distribution ${ \mathcal N }\left(\mu ,\hat{\sigma }\right);$
  • (b)  
    compute the value of the importance function for the current position g_i = n(θ_i); 
  • (c)  
    compute the target function f(θ) for the current position in the parameter space f_i = f(θ_i).
• 3.  
  The integration result is calculated as $\tfrac{1}{N}{\sum }_{i=1}^{N}\tfrac{{f}_{i}}{{g}_{i}}$.

Here, importance sampling is used to improve the convergence of the Monte-Carlo integration. The form of the specific importance function does not have any implication for the posterior PDF. Therefore, though we use the multivariate Gaussian as the importance function, the posterior PDF can still be an arbitrary function more or less confined in the parameter space.

One of the most frequent applications of the Bayesian analysis and MCMC is to infer parameters θ of a model M which is an analytical function that describes theoretical dependence of y upon x and has a set of free parameters θ:

$$\begin{eqnarray*}&&y=M(x,\theta )\end{eqnarray*}$$

from the observed data points (D = [X_i, Y_i]: i = 1..N) where N is the number of data points, and X_i and Y_i are empirically determined values of x and y in the ith measurement. The uncertainties of the fitted parameters $\theta =[{\theta }_{1},{\theta }_{2},\ \cdots ,\ {\theta }_{{N}_{p}}]$ also have to be estimated. SoBAT contains the (MCMC_FIT) routine which aims to solve this problem.

MCMC_FIT utilizes the assumption that the error corresponding to Y measurements is normally distributed with the standard deviation σ_Y. Thus, the likelihood function is the product of N Gaussians

$$\begin{eqnarray}&&P(D| \theta )=\displaystyle \frac{1}{{\left(2\pi {\sigma }_{Y}^{2}\right)}^{\tfrac{N}{2}}}\displaystyle \prod _{i=1}^{N}\exp \left\{-\displaystyle \frac{{\left[{Y}_{i}-M(X,\theta )\right]}^{2}}{2{\sigma }_{Y}^{2}}\right\}.\end{eqnarray} \tag{ 10 }$$

The measurement error σ_Y is considered as one of the unknown parameters. It is also assumed to be the same for all data points and is inferred during the MCMC simulations together with θ.

As a priori knowledge, a user can provide a range of the possible model parameter values θ:

$$\begin{eqnarray*}&&{\theta }_{i}^{\min }\leqslant {\theta }_{i}\leqslant {\theta }_{i}^{\max }.\end{eqnarray*}$$

Thus, our prior probability distribution can be expressed as follows

$$\begin{eqnarray}&&P(\theta )=\displaystyle \prod _{i=1}^{N}H({\theta }_{i},{\theta }_{i}^{\min },{\theta }_{i}^{\max }),\end{eqnarray} \tag{ 11 }$$

where $H({\theta }_{i},{\theta }_{i}^{\min },{\theta }_{i}^{\max })$ is the PDF of a uniform distribution in the range $[{\theta }_{i}^{\min },{\theta }_{i}^{\max }]$ which is defined as

$$\begin{eqnarray}H({\theta }_{i},{\theta }_{i}^{\min },{\theta }_{i}^{\max })=\left\{\begin{array}{cc}\tfrac{1}{{\theta }_{i}^{\max }-{\theta }_{i}^{\min }}, & {\theta }_{i}^{\min }\leqslant {\theta }_{i}\leqslant {\theta }_{i}^{\max }\\ 0, & \mathrm{otherwise}\end{array}\right..\end{eqnarray} \tag{ 12 }$$

One of the ways to check the correctness of the parameter inference is to estimate the posterior predictive distribution by sampling from it during the main sampling procedure. In the MCMC_FIT routine, Equation (10) is used to generate a sample from the posterior predictive distribution of the measured data [Y] for every sample from the posterior distribution [P(θ∣D)]. The generated samples are returned in the ppd_sample keyword parameter. An example of using MCMC_FIT routine for sampling posterior predictive distribution can be found in Listing 2.

In the case of a user-supplied posterior PDF, the user has to use the lower level MCMC_SAMPLES routine and is responsible for simulating samples from the posterior predictive distribution and returning them in the ppd_sample keyword within the user-supplied IDL function computing posterior PDF.

To test the sampling procedure used in the developed code, we selected the following 1D distributions: slightly asymmetrical triangular

$$\begin{eqnarray*}f(x)=\left\{\begin{array}{ll}\tfrac{2(x-a)}{(b-a)(c-a)} & \mathrm{for}\ a\lt x\leqslant c,\\ \tfrac{2(b-x)}{(b-a)(b-c)} & \mathrm{for}\ c\lt x\lt b,\\ 0 & \mathrm{otherwise}\end{array}\right.\end{eqnarray*}$$

with a = 0.5, b = 3, and c = 2.5 (see Figure 1(a)); uniform

$$\begin{eqnarray*}f(x)=\left\{\begin{array}{ll}\tfrac{1}{(b-a)} & \mathrm{for}\ a\lt x\lt b,\\ 0 & \mathrm{otherwise}\end{array}\right.\end{eqnarray*}$$

with a = 0.5 and b = 3 (see Figure 1(b)); exponential

$$\begin{eqnarray*}f(x)=\left\{\begin{array}{ll}\lambda {e}^{-\lambda x} & \mathrm{for}\ x\gt 0,\\ 0 & \mathrm{otherwise}\end{array}\right.\end{eqnarray*}$$

with λ = 1 (Figure 1(c)); and a bimodal mixture of two normal distributions with different expected values and dispersions

$$\begin{eqnarray*}&&f(x)=0.8\displaystyle \frac{1}{2\pi {\sigma }_{1}^{2}}{e}^{\tfrac{x-{\mu }_{1}}{2{\sigma }_{1}^{2}}}+0.2\displaystyle \frac{1}{2\pi {\sigma }_{2}^{2}}{e}^{\tfrac{x-{\mu }_{2}}{2{\sigma }_{2}^{2}}}\end{eqnarray*}$$

with μ₁ = 0, μ₂ = 7, σ₁ = 2, and σ₂ = 1 (see Figure 1(d)). Normalized histograms of the 10⁵ MCMC samples generated for each distribution are shown in Figure 1. The obtained histograms perfectly coincide with the corresponding target densities shown in Figure 1 with solid black lines.

Zoom In Zoom Out Reset image size

To demonstrate the correctness of the sampling procedure in the multiparametric case, we present the testing results for a set of bivariate target probability densities. We selected 2D versions of the distributions used in Section 4.1: pyramid (Figure 2(a)), 2D uniform distribution bounded by a square (Figure 2(a)), 2D exponential distribution, and a mixture of three bivariate normal distributions with different expected values and covariance matrices. The 2D histograms (see Figure 2) perfectly coincide with the target densities, shown in Figure 2 by contours.

Zoom In Zoom Out Reset image size

Let us consider a simple example of fitting a set of synthetic data points X_i, Y_i by a linear function to illustrate the practical usage of SoBAT. The synthetic data points in our example are generated using the linear dependence with the addition of normally distributed noise

$$\begin{eqnarray*}&&{Y}_{i}={{kX}}_{i}+b+N(0,\sigma ),\end{eqnarray*}$$

where k = 0.5, b = 1, and σ = 2.

First, we need to specify the model as a function describing the linear dependence of y upon x. The model function for the linear dependence is given in Listing 1.

Listing 1. 

Then, we define priors and initial guesses for the model parameters k and b (lines 2–6 in Listing 2). For the parameter k, the prior is defined as a normal distribution with zero expectation and standard deviation of 2, while for the parameter b, we set a uniform prior. After the call of MCMC_FIT function (lines 12–16 in Listing 2), the variable fit will contain the best-fitting values for Y. The fitted parameters values and corresponding uncertainties will be stored in the pars and credible_intervals variables. The MCMC samples will be returned in the samples keyword. The latter can be used to plot histograms approximating the marginalized posterior distributions. The histograms obtained for the slope (k), bias (b), and noise level (σ) are given in Figures 3(b)–(d). Note that the true parameter values (green vertical lines in Figure 3) do not coincide with global maximum of the histograms, but lie within the high probability area illustrated by the histograms. Such a behavior is expected because our inference (as any measurement) is uncertain. The uncertainty is described by the width of the histograms and can be quantified for an arbitrary level of significance by computing credible intervals as percentiles of the samples generated with the MCMC code.

Zoom In Zoom Out Reset image size

The keyword variable ppd_samples (line 16 in Listing 2) will contain samples from the predictive posterior distribution (PPD) in the form of an array with dimensions of [N_d, N_s], where N_d is the number of data points and N_s is the number of samples. Figure 4 demonstrates a 2D histogram of the PPD overplotted with the data points. The plot in Figure 4 characterizes our inference as being correct and consistent with the data, since the entire high probability area predicted by the PPD is covered with data points on one hand, and there are no data points located in the white colored areas where the PPD probability is nearly zero on the other hand.

Listing 2. 

Zoom In Zoom Out Reset image size

To illustrate quantitative comparison of different user-defined models, we use the same synthetic data set as in Section 5.1 with the linear dependence contaminated by white noise. Now we attempt to fit it with a second model with the quadratic dependence:

$$\begin{eqnarray}&&y={kx}+b+{{cx}}^{2}+N(0,\sigma ).\end{eqnarray} \tag{ 13 }$$

Listing 3 shows the IDL representation of this model.

Listing 3. 

The MCMC Bayesian inference is done for both models and then the models are compared by calculating the Bayes factor. Figures 5 and 6 show the MCMC inference results for the quadratic model given by Equation (13). Though the best fits and posterior predictive distributions (see Figures 4 and 6) are very similar, the histograms of marginal posterior distributions are found to be significantly broader in comparison with the linear case. This demonstrates that the additional quadratic term does not improve the fit. The χ² and reduced χ² metrics are almost the same for both models (see Table 1) and do not show any significant advantage of one model against the other.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

SoBAT includes the MCMC_EVIDENCE function which allows us to calculate Bayesian evidences and hence the Bayes factor for comparing the models as described in Listing 4, where samples_l and samples_q are the MCMC samples simulated using the linear and quadratic models, respectively. The computed Bayes factor (28.8) indicates very strong evidence in favor of the linear model. This result is expected since we generated the synthetic data using the linear dependence with the background normally distributed noise.

Listing 4. 

Coronal loops are frequently observed to perform large-amplitude, rapidly damped, transverse oscillations when perturbed by events such as flares and coronal mass ejections. Their rapid damping is explained by resonant absorption which causes a transfer of energy from the kink mode to the torsional Alfvén mode (e.g., see the recent review by De Moortel et al. 2016). Pascoe et al. (2013) proposed a method to infer the transverse density profile in the oscillating coronal loop using the shape of the damping profile of the kink oscillation (Pascoe et al. 2012, 2015, 2016a, 2019; Hood et al. 2013). The method was first applied in Pascoe et al. (2016b) using a Levenberg–Marquardt least-squares fit to the data using the IDL code MPFIT (Markwardt 2009). It was extended in Pascoe et al. (2017a) to include additional physical effects and also use Bayesian inference. Pascoe et al. (2017c) also included the presence of a large initial displacement of the loop equilibrium position. A benefit of the MCMC approach is that we can readily extend our models in this way, allowing us to investigate further details in the data.

We note that in previous applications of our MCMC code to coronal seismology (Pascoe et al. 2017a, 2017b, 2017c; Goddard et al. 2017), posterior summaries were given using the median value (and uncertainties by the 95% credible interval). Here (and in Pascoe et al. 2018, 2020) the MAP estimate is used rather than the median.

In this paper, we use the simplified version of the oscillation profile model published in Pascoe et al. (2017a):

$$\begin{eqnarray}y(t)={y}_{\mathrm{tr}}(t)+\left\{\begin{array}{ll}{A}_{0}{e}^{-{\left(\tfrac{\tilde{t}}{\tau }\right)}^{n}}\sin \left(\tfrac{2\pi \tilde{t}}{P}+{\phi }_{0}\right), & \tilde{t}\geqslant 0\\ {x}_{0}, & \tilde{t}\lt 0\end{array}\right.,\end{eqnarray} \tag{ 14 }$$

where ${\phi }_{0}=\arcsin \left(\tfrac{{x}_{0}}{{A}_{0}}\right)$ is the initial phase, A₀ is the initial amplitude, t₀ is the start time of the oscillation, $\tilde{t}=t-{t}_{0}$, P is the oscillation period, and x₀ is the initial displacement which prescribes the oscillation phase. The parameter n prescribes the damping profile. The background trend (y_tr) prescribes the equilibrium position and is calculated using spline interpolation from the reference points located at the time instances when the loop comes through the equilibrium (blue diamonds in Figure 7). The positions of the reference points are free parameters of the model and are identified during the Bayesian inference.

Zoom In Zoom Out Reset image size

As an example, we consider the time series of the loop position taken for Event 43 Loop 3 from the catalog of oscillations by Goddard et al. (2016). This loop is also referred to as Loop #1 in the seismological analysis by Pascoe et al. (2016b, 2017a). The observational data points and the best fit obtained using the MCMC_FIT function are shown in Figure 7. The histograms approximating marginal posterior distributions of oscillation period, amplitude, decay time, initial displacement, start time, and the position of a trend reference point are given in Figure 8.

Zoom In Zoom Out Reset image size

The posterior predictive distribution inferred using our MCMC code is given in Figure 9. The shaded area demonstrates the region on the plot where the data points are predicted to be observed. For a data consistent inversion, the measured data points should be located inside the shaded region and the shaded area itself should not broaden far away from the data points. That means that a model should predict the observed data points, but it should not predict observations being far away from the actually observed data.

Zoom In Zoom Out Reset image size