{'fraction_non_alphanumeric': 0.024225649956462248, 'fraction_numerical': 0.0032964299042169422, 'mean_word_length': 4.47166921898928, 'pattern_counts': {'":': 1, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bayesian parameter estimation of miss-specified models

Fitting a simplifying model with several parameters to real data of complex objects is a highly nontrivial task, but enables the possibility to get insights into the objects physics. Here, we present a method to infer the parameters of the model, the model error as well as the statistics of the model error. This method relies on the usage of many data sets in a simultaneous analysis in order to overcome the problems caused by the degeneracy between model parameters and model error. Errors in the modeling of the measurement instrument can be absorbed in the model error allowing for applications with complex instruments.


I. INTRODUCTION
A. General problem Whenever phenomena of nature are not directly accessible, it is a common way to build a model of them. These models usually have parameters for which one has to find the best fitting parameter values for each phenomena. Unfortunately, in this process of fitting and modeling there are two main problems. On the one hand, models are often derived from theoretical analysis of the underlying processes. However, these theoretical models describing complex real world systems are necessarily idealizations. The implied simplifications allow us to focus on the essential physics, to keep the model computationally feasible, and to investigate systems for which not all of their components are perfectly known. On the other hand there is always measurement noise in the observation process.
By judging theoretical models on a goodness-of-fit basis they often perform insufficiently when faced to real data of the actual system. The difference between real and modeled data can easily exceed the error budget of the measurement. The reason is that the idealized models do not capture all aspects of the real systems. But aspects, which are not modeled, imprint on the data as well. To discriminate these from measurement errors or noise, we use the term model errors to describe these imperfections of our theoretical description. A common approach to determine the model parameters from the data is a likelihood based methodology (X 2 -fit [1,2], maximum likelihood [3], or Bayesian parameter estimation [4]) that uses the measurement uncertainties as a metric in data space. The resulting parameter estimates can be strongly distorted by the desire of the method to minimize all apparent differences between predicted and real data, indifferently if these differences occur due to measurement or model errors.
Thus, the model error should be included into the error budget of the parameter estimation as well. However, this requires a model for the not yet captured aspects of the system or at least for the model errors which are produced by the missing aspects. In a rigorous way, we need to do a case by case analysis of the missing physics. However, this would be completely impractical in cases like the one we have in mind -spectroscopic investigation of complex astrophysical objects. Instead, this paper aims to construct a plausible, but by no means perfect, effective description of the model errors and their uncertainty. We want to show how the estimation of the model parameters is improved by taking a model error into account.
B. Example of the general problem in astrophysics: Spectroscopic investigation of complex astrophysical objects "Remember that all models are wrong, the practical question is how wrong do they have to be to not be useful." [5]. This quote of George E.P. Box, a frontier in Bayesian inference, is well applicable for astrophysical models. An example for the process of modeling and fitting is the task to fit stellar models to observed stellar spectra. Many different physical conditions determine the spectrum of a star e.g. stellar parameters like the effective temperature T eff , the magnitude of the gravitational surface acceleration log g, the metalicity [Fe/H] and abundances for almost the full periodic table. A usual approach to determine stellar conditions is fitting data to synthetic model spectra. The calculations of synthetic spectral models are normally done by radiative transfer calculations. However, if one accounts for all physical conditions in these calculations one can not keep the models computationally feasible. Thus, some simplifications have to be done e.g. local thermodynamical equilibrium (LTE) or plane parallel geometry. In the resulting star models only the most important physical parameters are included. We will denote physical conditions included in the process of modeling with p, while not modeled physical conditions will be denoted with c. Additionally, it is still very hard to produce high quality synthetic data due to inaccuracies in atomic and molecular constants. Thus, it is almost impossible to ensure spectral lines to fit perfectly over all stellar parameters and spectral ranges. Due to the inaccuracy in the modeling process stellar spectra can only be described approximately. As an illustrative, simple model for an absorption spectrum accounting for a gray body spectrum with absorption lines, one could use the following: where G(ν − ν i , σ 2 i ) denotes a Gaussian in ν with mean ν i and standard deviation σ i . The exponent of the powerlaw α, the effective temperature T eff and the strength parameters for the absorption line spectrum t 1 , t 2 , t 3 are assumed to be the unknown parameters, while the amplitude A, the mean frequency ν 0 , the means of the absorption peaks ν 1 , ν 2 , ν 3 , and the line widths σ 2 1 , σ 2 2 , σ 2 3 are known. This simplyfing model does not account for the gravitational surface acceleration log g, the metalicity [Fe/H] and all abundances of the periodic table expect for three. Assuming one determines the correct values p of the star models, not modeled conditions will present themselves in a model error. However, we have observations of many stars. If some of them are fitted with the same model, this gives us a chance to at least identify and characterize from star to star varying spectral structures due to model errors. Then, the model errors can be statistically added to the error budget. However, parameter determination and determination of the statistics of the model errors are coupled, which is the reason why they must be performed together. Here we use a generic, hierarchical Bayesian model based on Gaussian processes to describe the model errors, which is schematically shown in fig. 1. The structure of this work is the following. In sec. II we explain the general problem of parameter estimations of miss-specified models and give an overview of the state of the art methods tackling this problem. A brief discussion of the Wiener filter given in sec. III, since the understanding of it is crucial for our algorithms. In sec. IV we state and explain our model for the model error and its uncertainty and furthermore discuss our prior assumptions. In sec. V we derive and present the inference algorithm and show its capability on mock data. Finally, we compare the inferred results with a least-square parameter estimation.

II. PROBLEM SETTING AND STATE OF THE ART
Computer models simulating star spectra are invaluable tools in astrophysics. Unfortunately, there is typically uncertainty in the predictions of the synthetic star spectra models. The arising errors can be categorized in two main sources. The first one results from the use of an inadequate model with a correct set of parameters values. This means one has found the true parameter values for the wrong model, which is called the model error. The second one is the uncertainty in the parameters considering the true model, which is called the parameter error. By modeling and computing stellar spectra both errors are present simultaneously. Let us denote with t * the true forward operator of the model M and with t another forward operator, e.g. the forward operator used for the parameter estimations. Let p * be the true set of parameter values and p another set of parameter values. Then the parameter error is the difference between t * (p) and t * (p * ) which is independent of t, but a quantification would require the knowledge of t * . However, t * is not available to us as it represents the true physics. This raises the question of a meaningful parameter error evaluation, especially in the absence of a model error model. The direct evaluation of model errors, defined as is a subject of large uncertainty and has only little use in practice. The importance of model errors in parameter estimation is widely recognized, but there exists only little research about it. Thus, still no general methodology for quantifying the impact of model error for parameter estimations and its influence on parameter error analysis is available.
In literature a common way to describe the forward problem used to estimate the unknown parameters is described in the following way where d is the observed data, t(p) the forward operator belonging to the mathematical model M and e is the error. The error e is the sum of a measurement error n and a model error u e = n + u .
The inverse problem is formulated by finding the best fitting parameters p. Usually this is achieved by minimizing the residual e via a classical or weighted least-square method. In this method p is obtained by minimizing the sum of the squared differenceŝ where E represents the covariance of the error e. The underlying assumption for this formulation of least-square is that the measurement error and model error have the same statistics. However, there are two main problems.
On the one hand the model error is not random, but it is systematic and certainly will not have the same statistic as the measurement error. The model error will for different measurements of objects, which are described by the same model and same parameters, have a comparable strength on different positions and a certain correlation structure. On the other hand the covariance of the model error is not known a priori and can only be estimated. Thus, this approach is almost useless, if the model error u is dominating the overall error e.
Another main problem in estimating parameters and model error is that several approaches do not give a unique solution for the parameters p due to the degeneracy between the model error and the parameters. One way to tackle the problem of degeneracy is adding prior information on the parameters. The objective function minimized to find the best fitting parameters reads where P is the covariance matrix for the parameters and p e is a estimated set of parameter values. The assumption for the parameters is that they are Gaussian distributed. This approach allows to derive a posterior uncertainty of the parameters p, where all plausible sets of parameter values should be represented. Unfortunately, the model error is again not incorporated directly in the modeling process for the data observation and this approach has the same drawbacks as least-square.
Assuming we know the measurement error distribution, we can write the solutions for the inverse problem for d = t(p) + u + n as a set where s is the quantity of interest. This shows that the degeneracy of the model error and the parameters can be embedded on a manifold, on which the solutions are equally distributed without prior knowledge. The authors of [6] investigated the role of the priors on the model parameters and model errors in the inference. For one data set, all points of M s have equal likelihood and the posterior distribution on this manifold is completely determined by the prior. The probability to obtain the true values of p and respectively u will be high, if the prior has a high density for the correct solutions relative to the other solutions. In general one can distinguish two extreme cases.
For the first case, the authors of this study considered a weak prior on p and a diffuse, zero-mean prior distribution for u. It was found by [7], that the posterior will be flat over the set of (p, u) that result in the same s. As the authors of [7] assumed a zero-mean prior on u, the probability in regions, where u is near zero, will be higher. Thus, the estimated model parameters p will be similar to a maximum a posteriori estimator, which may not cover the true values of p.
The second case is given, if realistic prior information for u and/or p is available, which leads to a tighter posterior distribution. Unfortunately, if the prior information is wrong, a low probability is attached to the true value of p and the posterior uncertainty may fail to cover them. If a strong prior for the parameters is available from other observations that do not rely on the model process, the quality of model error and model parameter estimation is higher.
As we have seen, there are many subtleties in these methods. Thus, we will now discuss the main problems of approaches in literature to tackle the issue of a missspecified model and give a conceptual overview of our ideas. If we try to isolate the effect of the model errors for the parameter estimation with the help of eq. 7 this will lead to following problems.
The first problem is that the model error is not random and does certainly not have the same probabilistic statistics as the measurement error. To solve this issue we will try to model the model error as its own quantity with its own statistic.
The second problem is that although the inclusion of prior information on the parameters in the objective function reduces the ill-posedness of the inverse problem, it does not guarantee the uniqueness of the solution. As we are in the comfortable situation to have more than one data set with approximately the same parameter values for our analysis, the ill-posedness can be reduced further.
The third problem is that each data set can be fitted in different ways. Thus, one parameter set can be the best fit for the first data set and another set of parameter values can be the most appropriate for the next data set. We will analyze all available data sets at the same time, but do not assume that all data sets have the same true parameters. Note that this is not the same as combining the results of an analysis of each data set on its own, since the different data sets share the same model error statistics.

III. WIENER FILTER THEORY
In the following, we want to give a short overview of the Wiener filter. Stellar spectra, which we will also call the signal s in this work, are physical fields. Thus they have infinitely many degrees of freedom, which motivates the usage of the framework of information field theory (IFT) [8] in which the Wiener filter appears as the solution to the free field case. We observe a single stellar spectrum with the help of a measurement instrument that somehow discretizes a field, e.g. by integrating and averaging. In our mathematical description we represent the measurement instrument as a linear operator R. Further, we assume the instrument to be well known, such that errors arising from the modeling process of the measurement instrument will be ignored. The measured data d can be expressed by applying R to the signal s corrupted with noise n. The measurement equation describing this reads In particular, we assume the noise n to be Gaussian, i.e.
n G(n, N ) (10) where we denote with n † the complex conjugate and adjoint of n and N is representing the covariance of the noise distribution. Mathematical conventions in this formula, e.g. how an operator acts on a field and how the scalar product is defined, are explained in appendix A 1. By using the measurement eq. 9 the likelihood for the data given the signal and noise realization can be derived. The conditional probability P(d|s, n) = δ[d − (Rs + n)] tells us how probably the data is given the noise and the signal: Since we do not know the noise, we marginalize it out: For the task of infering the signal s from the data d, we have to inspect the posterior distribution. This is provided by Bayes theorem [9], P(s|d) = P(d|s)P(s) P(d) .
In analogy to statistical physics, Bayes theorem is reformulated by introducing the information Hamiltonian H(s, d) and the partition function Z(d) We will assume a Gaussian prior for the signal s, i.e.
where for the moment the prior covariance S of the signal s is assumed to be known. Together with the data model and Bayes equation, this allows us to calculate the information Hamiltonian H(s, d), which is given by A quadratic completion leads to a more conventional form of the information Hamiltonian, where we drop terms not depending on s, as denoted by the equality up to constants " =": As mentioned before, the posterior is the quantity of interest in the inference and with a short calculation we get Since the l.h.s and the r.h.s. of eq. 24 are normalized, holds exactly. Here we introduced the mean m = DR † N −1 d, which can be written as applying the infor- is known as the Wiener filter [10]. In many real world cases the prior uncertainty structure of the signal s and additional fields are not known a priori and one would like to obtain posterior uncertainties. One way to counteract these problems is proposed in [11,12] which reconstructs the uncertainty structure alongside the signal simultaneously. We will adopt this improved formulation in the following.

IV. MODEL ERROR MODEL
In this section we motivate our model for the model error and its uncertainty. Furthermore, we want to derive the information Hamiltonian for this specific choice of model error uncertainty. As our goal is to identify the model parameters of a theoretical model, we have to build a model for the model error. Additionally, we have to connect the model error to observed data. Assuming we have spectral observations of n stars, we denote the spectral data of star i with d i and the data vector d = (d i , · · · , d n ) describes the ordered collection of all data sets in the data vector. The measurement equation is given in eq. 9, where the signal can be described as a parameter model t(p) corrupted with a model error u. However, in the case of a spectroscopic analysis of an absorption spectrum, we want to ensure, that the signal is positive. Therefore, we adopt the following measurement model: This change to a positive definite signal has two main consequences. On the one hand the model error e u is now multiplicative to the model e t(p) , which is often a more accurate description than an additive model error as it allows a better separation of systematic model error and measurement error [13]. On the other hand, to compare the model to data, we have to take the logarithm of the model. In case of our spectroscopic model in eq. 1, this is given by The consequence is that now the model is linear in every parameter except for T eff . For a rigorous statistical analysis we need a prior P(p), which is also necessary to reduce the degeneracy between the parameters and the model error.
However, we also have to build a model for the statistics of the model errors u. We assume the model error to be Gaussian with a model error uncertainty covariance U . By analyzing data sets with the same model forward operator t and the same model-generating process, we can assume that the model errors u i follow the same statistics. For the model error uncertainty covariance U we try to be as general as possible for the reasons explained in sec. II. On the one hand, we have to take care of locally varying errors, e.g. peaks in stellar spectra, which are not modeled. This task can be tackled with a diagonal matrix e θ in position space (here ν, parameterized by the field θ). On the other hand we think that for some part of the model error budget a priori no position should be singled out. Thus, statistical homogenuity is assumed, which is according to the Wiener-Khintchin theorem [10,14], represented as diagonal matrix e τ in Fourier space (here the Fourier domain conjugate to ν). For signals with more than one dimensions statistical isotropy might be assumed as well. This means to assume all directions to be statistically equal, resulting only in a dependence of that part of the model error correlation on the relative distance of the compared locations. Thus, statistical isotropy leads to a collapse of the Fourier-diagonal part of the correlation structure to a one-dimensional power spectrum function e τ = e τ (κ) of the harmonic modes κ. This is incorporated in the uncertainty through the power projection operator P, which is explained in appendix A 2, and the Fourier transform operator F. One simple way to combine both parts of uncertainty structures into the uncertainty covariance is a modified product of both, which is given by The index y = ln(k) indicates that we assumed a logarithmic scale for the Fourier space coordinate of τ and the index x denotes a linear space coordinate for θ. We indicate the transformation of a field to a diagonal operator in the domain of the field by " ". In this parameterization the model error uncertainty is a multiplicative combination of a mask in position space e θx and a convolution. Any field processed by U (and also by its inverse) first goes through the (inverted) mask, then gets convolved (deconvolved) and at last the (inverted) mask again imprints on the field. Again for the sake of a statistical rigorous analysis, we need a prior for τ and θ. To account for the smoothness of the power spectrum on a logarithmic scale, a smoothness prior is adopted to τ . Similarly, a smoothness prior for θ on a linear scale in position space is adapted to θ: The strength parameters of the curvature, σ θ and σ τ , represent the tolerated deviation by weighting the L 2norm of the second derivative. For the limit σ θ\τ → ∞ the expected roughness is very high and therefore rough realizations are not suppressed. Otherwise, for the limit σ θ\τ → 0 the smoothness prior enforces a linear behavior in the chosen coordinate systems, which correspond to a power-law for the power spectrum. However, the exact values of σ θ\τ are not the driving forces for the inference.
Since the priors on θ and τ only punish curvature, there is a degeneracy in U . We can add a constant factor a ∈ R to θ, if we subtract at the same time 2a from τ . This operation leaves the functional form of U invariant. As one is not interested in the individual values of τ or θ, this does not affect the Hamiltonian. One can use the degeneracy to control the numerical stability of the algorithm by keeping θ and τ in numerically stable regions.
In the following the model error parameters (τ, θ) will be merged into the vector q = (τ, θ). Thus, if we talk about parameters without further specification, both p as well as q are meant. The hierarchical Bayesian network summarizing our assumptions can be seen in fig. 2.
Similar as for the Wiener filter theory in sec. III, we introduce the posterior information Hamiltonian H(p, τ, θ, u|d) and the partition function Z(d): The resulting full Hamiltonian can be calculated as before, such that we get (36) The model error u is assumed to be generated by drawing from a Gaussian distribution, which is neither diagonal in position nor in Fourier space. Furthermore, it is possible to draw from a Gaussian distribution in its eigenbasis by weighting an independent, white spectrum excitation field ξ with unit variance with the square root of its eigenvalues. We will adopt this idea and reformulate our problem in the following way.
Here we introduced the amplitude operator A. The change in causality can be seen in fig. 3, where one observes that now all unknown fields directly interact with the data. This reparametrization is following the ideas of [15,16]. Inserting the change of variables from u to ξ into the Hamiltonian gives   IV. The parameters in the top row, στ and σ θ , as well as P(p) have to be specified by the user. The parameters characterizing the model error uncertainty θ, τ and the model parameters p will be infered by the algorithm.
All terms that are unimportant for the later process are collected in H 0 : (39) These terms will not affect the further inference. For the inference we will need the posterior Hamiltonian H(p, τ, θ, ξ|d) = H(p, τ, θ, ξ, d) with ln [Z(d)] defined in eq. 35. Since the partition sum does not depend on ξ, q and p, we can ignore it in the following. In this section we want to present the approximation of the posterior distribution via variational Bayes. Thus, we adapt the methods presented in [11,12] which performs a minimization of all unknown fields at the same time and provides a posterior uncertainty estimation via samples.

A. Approximation of the posterior distribution
To use the implementation in the Numerical Information field theory package (NIFTy) [17,18], we have to ensure that all used priors, except the noise prior, are white priors. Due to this fact, we have to reformulate our parameterization of U in eq. 30. As already mentioned the reformulation is possible by weighting a white excitation field with the square root of the covariance. This can be done for the smoothness priors on θ and τ .
The exact transformations and details on the chosen coordinates for the excitation fields are given in appendix A 3. Thus, the covariance operator U has the reformulated form Now, we can again use a variable transformation to write The coordinate transformation of the parameters p into a coordinate system ξ p where the prior is white is a little bit more subtle, since the prior for the parameters can almost have any form. However, as long as the cumulative distribution function of the prior can be calculated analytically, this coordinate transformation can always be applied [15,16].
As a next step one stores all unknown fields, namely ξ, ξ θ , ξ τ , ξ p , in a multi-field φ. Now, we can continue as before and calculate the full Hamiltonian (44) Usually one would need to calculate the derivatives of this in order to be able to apply a Newton scheme. However, since the model for the model error is constructed from NIFTy model objects, NIFTy is able to compute all derivatives regarding ξ, ξ θ , ξ τ on its own. Only the gradient with respect to the model parameters ξ p has to be calculated and implemented, which is individual for each model t(ξ p ). In order to apply a variational inference method, the posterior distribution has to be approximated. The advantage of the method implemented in NIFTy is that we can approximate the posterior as a Gaussian in all unknown fields whereφ represents the mean of the multi-field φ and Φ is its full covariance matrix, which contains the crosscorrelations of the different fields within the multi-field. By using this approximation one can minimize the variational Kullback-Leibler divergence (D KL ), which is given by

B. Numerical example
We demonstrate the applicability of our derived algorithm for one-dimensional synthetic data. The algorithm is implemented in python using the NIFTy package. By comparing our algorithm with the least-square algorithm, we show its capability. In our example the signal has a resolution of 128 pixels. As a prior for the parameters we used a Gaussian, such that where we assumed the prior matrix P −1 to be diagonal. As discussed before, we have to transform the parameters, such that they have a white prior. For the chosen prior this is straight forward and the transformation is given by where p * i is the mean estimation of the i-th parameter and σ i = P . Thus, we can draw a excitation field ξ p from the white prior. The corresponding models which are used in the example with 32 different data sets and models can be seen in fig.  4. We produce the data sets according to eq. 27, where the signal field s is a sum of the error and the model error. By drawing a white excitation field ξ θ and fixing the power-spectrum to we construct the operator A. The noise covariance is given by 0.2×1 and the measurement instrument is a unit response. The gradient of the model with respect to the transformed parameters Fig. 5 shows one of the spectra with the corresponding model. The model with its free parameters is not able to capture all the aspects of the spectrum. Conventional algorithms would try to explain the differences through the model, whether they are due to the model error or to the model. In fig. 6 the 32 different model errors and the square root of the covariance U are shown, which illustrates the correlation structure and the difficulty to reconstruct a spatial varying uncertainty covariance. One can observe that for this special U the model error mainly imprints on the spectra in their ascending parts and the part where the absorption lines are present. Fig. 7 shows a comparison on between our method and a least-square method (LS) for the estimation of model parameters. Using multiple records can cause the optimization scheme to go to a local minimum. After some tests, we decided to divide the likelihood Hamiltonian by the number of data sets used in the particular inference similar to the proposed solution of this problem by [19].. We can observe that our method is able to learn from an increasing number of data sets while the LS estimates each data set on its own and thus does not improve with more data sets. However, our method saturates at about 512 data sets on an average parameter error of 6% and the addition of more data sets does not further reduce the error. As for each model still the parameters p i have to be found individually in the presence of the unknown noise relations n i and u i . The square root of the diagonal of the model error uncertainty covariance U and its reconstruction together with the root mean square error u RMS are shown in fig. 8. One would expect that the recovered uncertainty, as an conservative estimator, is in general smaller than the correct uncertainty. However, in our reconstruction of the model error uncertainty the peaks can be modeled well, while the tails are overestimated.

VI. CONCLUSION
Estimations of parameters from models are essential to get insights into physics. However, many physical objects are not understood in their full complexity. Thus, the model errors need to be included in the error uncertainty, but this would afford a case by case analysis of the missing physics and therefore is often impractical. To address these problems this paper constructs a plausible and effective description of model errors and their uncertainty. To constraint the additional degrees of freedom it relies on using several data sets at the same time in a joint inference. We derive an algorithm to reconstruct the parameters of a model as well as the model error and its statistics from noisy data. By using more than one data set of the same type of stars, we got a more sophisticated estimate for the model error uncertainty. It was important to build a very general approximating description for the model error uncertainty in order to be capable for most sorts of model errors. One can see that a few general assumptions about the nature of model errors, i.e. locality and smoothness, can be sufficient to reconstruct the model error uncertainty when one has access to a large set of data. The algorithm can be graphically visualized through a Bayesian hierarchical model. It was able to reconstruct the model errors as well as the model parameters. Additionally, expectation means can be calculated through the samples from the inverse metric. By using approximate posterior samples the opportunity to deal with more complex uncertainty structures is computationally feasible. We can conclude that this paper enables the opportunity to deal with parameter estimates the spectral smoothness. The subscript lin indicates that the power space has linear bin bounds. We will denote the exponential of θ as A reformulation for the convolutional part of the model error uncertainty U already exists in NIFTy. This refor-mulation splits the τ field into two different fields, one describing the slope ξ φ and the other one describing the deviation of a linear slope ξ τ . This model works with white priors and thus is used here.