Introduction

To analyze a dynamic system both numerical simulations and experimental measurements can be performed. However, results of both approaches can deviate significantly. A possible reason for this mismatch might be the uncertainty in the input parameters for the numerical model (e.g. Young modulus, mass density). Therefore, a method to update the model parameters can be used so that the discrepancy between experiments and simulation is decreased. In this project Bayesian updating method is investigated.

Theory

In the following section some theoretical aspects related to Bayesian updating procedure are addressed.

Dynamics of the System

As described in Modal Analysis, free vibration of a linear elastic system is governed by the following equation:

$\begin{array}{l}\displaystyle \left( - \omega_i^2 \mathbf{M} + \mathbf{K} \right) \pmb{\phi}_i = \mathbf{0},\end{array}$

where $\begin{array}{l}\mathbf{M}\end{array}$ is the mass matrix, $\begin{array}{l}\mathbf{K}\end{array}$ is the stiffness matrix and $\begin{array}{l}\omega_i\end{array}$ and $\begin{array}{l}\pmb{\phi}_i\end{array}$ are the ith circular eigenfrequency and corresponding eigenvector. The system matrices can for example be obtained from the FEM discretization.

Now assume an experiment was performed and measured modal data were obtained. The dynamic equation can then be modified as follows:

(1)	$\begin{array}{l}\displaystyle \left( - \omega_{m,i}^2 \mathbf{M} + \mathbf{K} \right) \pmb{\phi}_{m,i} = \pmb{\varepsilon}_i,\end{array}$

where $\begin{array}{l}\omega_{m,i}\end{array}$ and $\begin{array}{l}\pmb{\phi}_{m,i}\end{array}$ are the ith measured eigenfrequency and ith measured eigenvector and $\begin{array}{l}\pmb{\varepsilon}_i\end{array}$ is ith error vector. Note that $\begin{array}{l}\pmb{\varepsilon}_i\end{array}$ is a vector, considering all measured eigenfrequencies and eigenmodes, eq. (1) results a matrix which components are $\begin{array}{l}\varepsilon_{ij}\end{array}$ . The latter matrix can be used as a measure of the mismatch between the measured data and modal properties computed based on the FEM discretization. In the case when $\begin{array}{l}\omega_{m,i}\end{array}$ and $\begin{array}{l}\pmb{\phi}_{m,i}\end{array}$ coincide with $\begin{array}{l}\omega_i\end{array}$ and $\begin{array}{l}\pmb{\phi}_i\end{array}$ , which implies that numerical results correspond to the measurements, components $\begin{array}{l}\varepsilon_{ij}\end{array}$ will be 0.

Bayesian Updating

The Bayesian updating method is based on the Bayes' theorem. It can be derived by applying the rule for conditional probabilities and is written as follows:

(2)	$\begin{array}{l}\displaystyle P(\mathbf{E}\|\mathbf{D})=\frac{ P(\mathbf{E}) \cdot P(\mathbf{D}\|\mathbf{E})}{P(\mathbf{D})},\end{array}$

where $\begin{array}{l}\mathbf{E}\end{array}$ is a vector of the parameters to be updated and $\begin{array}{l}\mathbf{D}\end{array}$ stands for a matrix containing the measured eigenfrequencies and eigenvectors. The term $\begin{array}{l}P(\mathbf{E})\end{array}$ is the probability distribution function (PDF) of the parameter vector $\begin{array}{l}\mathbf{E}\end{array}$ , which is also often called the prior distribution. $\begin{array}{l}P(\mathbf{D}|\mathbf{E})\end{array}$ means the probability distribution of the data for a given parameter vector, also denoted as a likelihood. $\begin{array}{l}P(\mathbf{E}|\mathbf{D})\end{array}$ is the posterior probability distribution function and expresses the updated probability distribution function of the parameter vector $\begin{array}{l}\mathbf{E}\end{array}$ , given the measured data. $\begin{array}{l}P(\mathbf{D})\end{array}$ is the normalization factor. In the subsequent sections some of the latter terms will be explained more detailed.

Prior Probability Distribution Function

Prior distribution function represents a prior belief about the probability distribution of the vector of parameters $\begin{array}{l}\mathbf{E}\end{array}$ in absence of any data. In this project a multivariate normal distribution was considered, the probability density function (PDF) is defined as follows:

(3)	$\begin{array}{l}\displaystyle f_{\mathbf{E}}(\mathbf{E})=\frac{1}{ \sqrt{\left\|\mathbf{\Sigma} \right\| \left( 2 \pi \right)^{k} }}} \mathrm{exp} \left[ {-\frac{1}{2}(\mathbf{E}-\pmb{\mu})^T\mathbf{\Sigma^{-1}(\mathbf{E}-\pmb{\mu})} \right],\end{array}$

where $\begin{array}{l}k\end{array}$ is the dimension of the vector of parameters, $\begin{array}{l}\pmb{\mu}\end{array}$ is the mean vector and $\begin{array}{l}\mathbf{\Sigma}\end{array}$ is the covariance matrix.

For further investigated examples an assumption of independent components of parameter vector $\begin{array}{l}\mathbf{E}\end{array}$ was taken into account. In this case, eq. (3) yields

(4)	$\begin{array}{l}\displaystyle P(\mathbf{E}) = f_{\mathbf{E}}(\mathbf{E}) = \frac{1}{Z_E(\alpha)} \mathrm{exp} \left[ -\frac{1}{2} \[\sum_{i=1}^{k} \alpha_i (E_i - \mu_i)^2 \] \right].\end{array}$

In eq. (4) coefficients $\begin{array}{l}\alpha_i=1/\mathrm{Var}(E_i)\end{array}$ were used to represent an inverse of the variance of the individual parameters and the normalization factor was abbreviated by $\begin{array}{l}Z_E(\alpha)\end{array}$ . Coefficients $\begin{array}{l}\alpha_i\end{array}$ are used to follow the notation of [1].

A drawback of the multivariate normal distribution is that also negative values of parameters $\begin{array}{l}\mathbf{E}\end{array}$ are allowed. If the individual components are material properties (e.g. Youngs modulus or mass density), negative parameters would yield nonphysical values. For such cases different prior distribution might be chosen. A convenient option is a lognormal distribution, which marginal PDF is defined as follows [2]:

(5)	$\begin{array}{l}\displaystyle f_E(E)= \frac{1}{E \sigma_{lnE}\sqrt{2\pi}}\mathrm{exp} \left[ -\frac{1}{2} \left( \frac{\mathrm{ln}E-\mu_{lnE}}{\sigma_{lnE}} \right) ^2 \right].\end{array}$

Here $\begin{array}{l}\mu_{lnE}\end{array}$ and $\begin{array}{l}\sigma_{lnE}\end{array}$ are the parameters of the lognormal distribution. Note that these parameters are mean and standard deviation of the natural logarithm of the variable $\begin{array}{l}E\end{array}$ . They can be computed from mean $\begin{array}{l}\mu_{E}\end{array}$ and standard deviation $\begin{array}{l}\sigma_{E}\end{array}$ of the variable $\begin{array}{l}E\end{array}$ itself based on the following formulas:

$\begin{array}{l}\displaystyle \mu_{lnE} = \mathrm{ln} \left( \frac{(\mu_E)^2}{\sqrt{(\mu_E)^2+(\sigma_E)^2}} \right) ,\ \sigma_{lnE} = \sqrt{ \mathrm{ln} \left( 1 + \frac{(\sigma_E)^2}{(\mu_E)^2} \right)}.\end{array}$

As the components of the vector $\begin{array}{l}\mathbf{E}\end{array}$ are assumed to be independent, the joint PDF is computed as a product of the marginal PDFs

(6)	$\begin{array}{l}\displaystyle f_{\pmb{E}}(\pmb{E}) = \prod_{i=1}^k f_{E_i}(E_i).\end{array}$

Likelihood Distribution Function

The likelihood term $\begin{array}{l}P(\mathbf{D}|\mathbf{E})\end{array}$ describes how well does a given parameter vector represent the measured data. In the scope of this project the likelihood function is defined with the use of the error matrix $\begin{array}{l}\varepsilon_{ij}\end{array}$ which was discussed in section Dynamics of the System. It reads:

(7)

$\begin{array}{l}\displaystyle P(\mathbf{D}|\mathbf{E}) = \frac{1}{Z_D(\beta)} \mathrm{exp} \left[ - \beta \[\sum_{i=1}^{k} \sum_{j=1}^{k} \varepsilon_{ij}^2 \] \right] = \frac{1}{Z_D(\beta)} \mathrm{exp} \left[ - \beta \[\sum_{i=1}^{k} \sum_{j=1}^{k} {\left[\left( - \omega_{m,i}^2 \mathbf{M} + \mathbf{K} \right) \pmb{\phi}_{m,i} \right]_j}^2 \] \right] ,\end{array}$

where $\begin{array}{l}\beta\end{array}$ is a coefficient expressing the contribution of the measured modal data and $\begin{array}{l}Z_D(\beta)\end{array}$ is the normalization factor which can be calculated as follows:

$\begin{array}{l}\displaystyle Z_D(\beta) = \int \mathrm{exp} \left[ - \beta \[\sum_{i=1}^{k} \sum_{j=1}^{k} {\left[\left( - \omega_{m,i}^2 \mathbf{M} + \mathbf{K} \right) \pmb{\phi}_{m,i} \right]_j}^2 \] \right] \mathrm{d} \mathbf{D}.\end{array}$

Note that in eq. (7) the components of the system matrices $\begin{array}{l}\mathbf{K}\end{array}$ and $\begin{array}{l}\mathbf{M}\end{array}$ depend on the parameter vector $\begin{array}{l}\mathbf{E}\end{array}$ , therefore the likelihood is a function of the parameter vector.

If the system matrices correspond to the measured modal properties, the components of the error matrix $\begin{array}{l}\varepsilon_{ij}\end{array}$ will be zero and for such parameter vector the likelihood term will be maximal.

Normalization Factor

Normalization term ensures that posterior function is a probability distribution and fulfills the normalization rule:

$\begin{array}{l}P(\mathbf{D}) = \int { P(\mathbf{E}) \cdot P(\mathbf{D}|\mathbf{E}) }\mathrm{d} \mathbf{E} = \int { \frac{1}{Z_E(\alpha)} \mathrm{exp} \left[ { -\frac{1}{2} \[\sum_{i=1}^{k} \alpha_i (E_i - \mu_i)^2 \]} \right] \cdot \frac{1}{Z_D(\beta)} \mathrm{exp} \left[{ - \beta \[\sum_{i=1}^{k} \sum_{j=1}^{k} {\left[\left( - \omega_{m,i}^2 \mathbf{M} + \mathbf{K} \right) \pmb{\phi}_{m,i} \right]_j}^2 \]}} \right] \mathrm{d} \mathbf{E} = \\ = \frac{1}{Z_E(\alpha) \cdot Z_D(\beta)} \int { \mathrm{exp} \left[{ -\frac{1}{2} \[\sum_{i=1}^{k} \alpha_i (E_i - \mu_i)^2 \] - \beta \[\sum_{i=1}^{k} \sum_{j=1}^{k} {\left[\left( - \omega_{m,i}^2 \mathbf{M} + \mathbf{K} \right) \pmb{\phi}_{m,i} \right]_j}^2 \]} }\right]\mathrm{d} \mathbf{E}\end{array}$

Posterior Probability Distribution Function

By substitution of the individual terms into the eq. (2) the posterior distribution function is obtained

(8)

$\begin{array}{l}P(\mathbf{E}|\mathbf{D})=\frac{ P(\mathbf{E}) \cdot P(\mathbf{D}|\mathbf{E})}{P(\mathbf{D})} = \dfrac{ \frac{1}{Z_E(\alpha)} \mathrm{exp} \left[ { -\frac{1}{2} \[\sum_{i=1}^{k} \alpha_i (E_i - \mu_i)^2 \]} \right] \cdot \frac{1}{Z_D(\beta)} \mathrm{exp} \left[ { - \beta \[\sum_{i=1}^{k} \sum_{j=1}^{k} {\left[\left( - \omega_{m,i}^2 \mathbf{M} + \mathbf{K} \right) \pmb{\phi}_{m,i} \right]_j}^2 \]} \right] } { \int \frac{1}{Z_E(\alpha) \cdot Z_D(\beta)} { \mathrm{exp} \left[ { -\frac{1}{2} \[\sum_{i=1}^{k} \alpha_i (E_i - \mu_i)^2 \] - \beta \[\sum_{i=1}^{k} \sum_{j=1}^{k} {\left[\left( - \omega_{m,i}^2 \mathbf{M} + \mathbf{K} \right) \pmb{\phi}_{m,i} \right]_j}^2 \]}} \right] \mathrm{d} \mathbf{E}} = \\ = \dfrac{ \mathrm{exp} \left[ { -\frac{1}{2} \[\sum_{i=1}^{k} \alpha_i (E_i - \mu_i)^2 \] - \beta \[\sum_{i=1}^{k} \sum_{j=1}^{k} {\left[\left( - \omega_{m,i}^2 \mathbf{M} + \mathbf{K} \right) \pmb{\phi}_{m,i} \right]_j}^2 \]} \right] }{\int { \mathrm{exp} \left[ { -\frac{1}{2} \[\sum_{i=1}^{k} \alpha_i (E_i - \mu_i)^2 \] - \beta \[\sum_{i=1}^{k} \sum_{j=1}^{k} {\left[\left( - \omega_{m,i}^2 \mathbf{M} + \mathbf{K} \right) \pmb{\phi}_{m,i} \right]_j}^2 \]} } \right] \mathrm{d} \mathbf{E}}.\end{array}$

If the number of parameters to be updated is large, the normalization integral in the latter equation might be difficult to compute. Therefore, numerical sampling method from the posterior distribution needs to be used instead.

Markov Chain Monte Carlo (MCMC)

The application of the previous Bayesian approach can be useful to update a model to predict a set of measured properties $\begin{array}{l}\bf{y}\end{array}$ from a system, e.g. modal properties, from a set of updated parameters $\begin{array}{l}\bf{E}_i\end{array}$ . A Monte Carlo technique generates as many samples of these vectors as needed, making it possible to get distributions of the predicted properties from a system. This is especially useful for estimations of the precision of a model with respect to a set of measurements or the evaluation of the error of that model. The objective is to get a set of updated parameters that will predict correctly the measured properties. The probability of $\begin{array}{l}\bf{y}\end{array}$ given $\begin{array}{l}\bf{D}\end{array}$ can be written as

(9)	$\begin{array}{l}\displaystyle P(\mathbf{y}\|\mathbf{D})=\int P(\mathbf{y}\|\mathbf{E})P(\mathbf{E}\|\mathbf{D})d\mathbf{E}}\end{array}$

Again, this equation is difficult to compute for a large number of parameters, as it depends on Equation (8), so a Markov Chain Monte Carlo (MCMC) approach can be used to approximate it. This iterative method is especially useful for the estimation of probability distributions, as it allows to obtain a distribution of the proposed parameters and its predictions.

Markov chain approach

MCMC is an iterative method based on the consecutive sampling from an initial proposed distribution from the parameters. The new sampled models are then accepted or rejected based on some specific criteria. This succession of models (or states) composes the so-called random walk of the method. If the chain is well constructed, the random walk will converge around an equilibrium distribution. Then the state of the chain after reaching the equilibrium is taken as the distribution of the model.

Illustration of Markov Chain Monte Carlo from [4]

In our case, the MCMC algorithm would proceed as follow:

A measured set of modal properties (eigenfrequencies and eigenmode shapes) is obtained. This can be from any source: experimental measurements, simulations, etc.
We set up a theoretical distribution of the $\begin{array}{l}\bf{E}\end{array}$ parameters to draw samples from.
A first state is calculated from a sample drawn from $\begin{array}{l}\bf{E}\end{array}$ .
A criterion for the acceptance of this set of parameters is applied. This criterion must ensure the convergence of the MCMC. In our case, we will use the Metropolis-Hastings algorithm.
If the sample is accepted, the distribution is updated with it.
Steps 3 to 5 are repeated for a set number of iterations or until some convergence criteria is met.

When the algorithm is finished, we can take the random walk as the distribution of the set of parameter. This distributions can be constructed based on those of the model parameters $\begin{array}{l}\bf{E}_i\end{array}$ or based on the modal properties $\begin{array}{l}\bf{y}_i\end{array}$ . Statistical analysis of these distributions can be later performed. The Equation (9) then approximates the vector of modal properties

(10)	$\begin{array}{l}\displaystyle \mathbf{\tilde{y}}\approx \dfrac{1}{L}\sum_{i=1}^{R+L-1}G(\mathbf{E}_i)\end{array}$

where $\begin{array}{l}G(\bf{E}_i)\end{array}$ is the outcome or the model, $\begin{array}{l}L\end{array}$ is the number of retained states and $\begin{array}{l}R\end{array}$ the discarded initial states ("burn-out" steps, see section Considerations on the convergence of the algorithm). Therefore, the mean of the distribution of the sampled modal properties is taken as an approximation of the real ones. If the method is constructed correctly this approximation will converge to the measured ones. However, the initial parameters $\begin{array}{l}\mathbf{E}_i\end{array}$ would not necessarily converge due to the contribution of the prior probability.

Metropolis-Hastings algorithm

The Metropolis-Hastings algorithm, or simply Metropolis algorithm, is one of the most used for generating a convergent Markov chain. It is based on the calculation of the probabilities expressed in the previous section. The algorithm is as follows:

(11)

$\begin{array}{l}If\text{ }P_{new}(\mathbf{E}|\mathbf{D})>P_{old}(\mathbf{E}|\mathbf{D}),\text{ } accept \text{ }state \text{ }\mathbf{E}_{new}\\ else,\text{ } accept \text{ }\mathbf{E}_{new}\text{ }with\text{ }probability\text{ }\dfrac{P_{new}(\mathbf{E}|\mathbf{D})}{P_{old}(\mathbf{E}|\mathbf{D})}\end{array}$

The use of the Metropolis algorithm ensures the convergence of the chain, although the speed of convergence is not known beforehand. However, on the contrary to some rejection sampling methods, the Metropolis algorithm doesn't suffer from the "curse of dimensionality". This means, that the rejection probability does not increase exponentially with the number of parameters, making it suitable for stochastic processes where a big number of parameters are involved. Some factors for the speed of convergence of the chain have to still be considered.

In some cases, the probabilities result in very small values that can be problematic to deal with floating points in a computer. For those situations, given that the probabilities for the chosen approach are usually modelled as exponential functions, the Metropolis algorithm can adopt a logarithmic form:

(12)

$\begin{array}{l}If\text{ } ln(P_{new}(\mathbf{E}|\mathbf{D}))< ln(P_{old}(\mathbf{E}|\mathbf{D})),\text{ } accept \text{ }state \text{ }\mathbf{E}_{new}\\ else,\text{ } accept \text{ }\mathbf{E}_{new}\text{ }with\text{ }probability\text{ }ln(P_{new}(\mathbf{E}|\mathbf{D}))-ln({P_{old}(\mathbf{E}|\mathbf{D})})\end{array}$

In this formulation, only the exponents of the probabilities have to be computed, being more suitable for high precision and small value computations. Also, sums are much less expensive than divisions from a computational point of view, which is relevant for iterative processes like MCMC.

Considerations on the convergence of the chain

First, we can observe that the first samples of the chain are quite unstable and usually diverge greatly from the objective equilibrium. If these samples are taken into account, the resultant distribution may have some outliers derived from those steps, which deteriorates artificially the quality of the distribution and delays reaching the equilibrium. In many cases, a first arbitrary number of accepted samples is disregarded, the so-called "burn-out" period. The final distribution is only taken into account after some stabilization of the chain is observed. This is especially relevant when the new model takes into account the standard deviation of the samples for doing the updating. In the case of this project, the standard deviation of the sampled model's distributions is kept constant, and the only update happens in the mean value of the distribution from where the parameters are sampled. Therefore, the benefits of having a burn-out period are negligible, so it will not be taken into account.

Another big factor in the convergence of the algorithm is the rejection rate and its relationship with the variance $\begin{array}{l}\sigma^2\end{array}$ of the distribution of the parameters. The rejection rate is the number of rejected samples with respect to the total number of proposed ones. According to Robert, Gelman and Wilk [5], the optimal rejection rate for parameters modelled as independent random normal distributions is 50%, being around 24% for co-joint normal distributions. A variance too big will result in a very high rejection rate, slowing the convergence due to the variability. If the variance is too low, the chain will stabilize very soon on a solution that is not the equilibrium one, so the convergence of the MCMC will be slowed. In case that the variances are prescribed by a real model, the orders of magnitude of the $\begin{array}{l}\alpha\end{array}$ and the $\begin{array}{l}\beta\end{array}$ have to be carefully chosen, not only to reflect the influence of their respective components but also for accelerating the convergence.

Pseudocode

For a better understanding of the procedure, a workflow of the code is described in this section.

Define the mean vector and covariance matrix of the prior distribution of the vector of parameters $\begin{array}{l}\mathbf{E}\end{array}$
Obtain measured eigenfrequencies $\begin{array}{l}\omega_{m,i}\end{array}$ and eigenvectors of the system $\begin{array}{l}\pmb{\phi}_{m,i}\end{array}$
Draw a sample $\begin{array}{l}\mathbf{E}_p\end{array}$ from the prior distribution
$\begin{array}{l}\mathbf{for} \ i = 1: \mathrm{N_{samples}}\end{array}$
1. $\begin{array}{l}\mathbf{if} \ i > 1\end{array}$
  1. Draw a sample from a proposal distribution $\begin{array}{l}\mathbf{E}_p \propto N(\mathbf{E}_a, \mathbf{I})\end{array}$
2. Compute system matrices $\begin{array}{l}\mathbf{K}\end{array}$ and $\begin{array}{l}\mathbf{M}\end{array}$ based on the sample $\begin{array}{l}\mathbf{E}_p\end{array}$
3. Evaluate the error matrix components $\begin{array}{l}\varepsilon_{ij}\end{array}$ from eq. (1)
4. Compute the prior $\begin{array}{l}P(\mathbf{E}_p)\end{array}$ and the likelihood $\begin{array}{l}P(\mathbf{D}|\mathbf{E}_p)\end{array}$ for the current sample $\begin{array}{l}\mathbf{E}_p\end{array}$ based on eq. (4) and eq. (7)
5. Evaluate the posterior (neglecting the normalization term) as: $\begin{array}{l}P_{new}(\mathbf{E}_p|\mathbf{D}) = P(\mathbf{E}_p) \cdot P(\mathbf{D}|\mathbf{E}_p)\end{array}$
6. $\begin{array}{l}\mathbf{if} \ i = 1\end{array}$
  1. Rename the sample: $\begin{array}{l}\mathbf{E}_a=\mathbf{E}_p\end{array}$
  2. Save the posterior value as: $\begin{array}{l}P_{old}(\mathbf{E}_p|\mathbf{D}) = P_{new}(\mathbf{E}_p|\mathbf{D})\end{array}$
7. $\begin{array}{l}\mathbf{else} \ \rightarrow\end{array}$ Apply Metropolis-Hastings algorithm
  1. $\begin{array}{l}\mathbf{if} \ P_{new}(\mathbf{E}_p|\mathbf{D}) > P_{old}(\mathbf{E}_p|\mathbf{D}) \ \rightarrow\end{array}$ sample accepted
    1. Rename the sample: $\begin{array}{l}\mathbf{E}_a=\mathbf{E}_p\end{array}$
    2. Save the posterior value as: $\begin{array}{l}P_{old}(\mathbf{E}_p|\mathbf{D}) = P_{new}(\mathbf{E}_p|\mathbf{D})\end{array}$
  2. $\begin{array}{l}\mathbf{elseif} \ P_{new}(\mathbf{E}_p|\mathbf{D}) / P_{old}(\mathbf{E}_p|\mathbf{D}) > \mathrm{rand}(0,1) \ \rightarrow\end{array}$ sample accepted
    1. Rename the sample: $\begin{array}{l}\mathbf{E}_a=\mathbf{E}_p\end{array}$
    2. Save the posterior value as: $\begin{array}{l}P_{old}(\mathbf{E}_p|\mathbf{D}) = P_{new}(\mathbf{E}_p|\mathbf{D})\end{array}$
  3. $\begin{array}{l}\mathbf{else} \ \rightarrow\end{array}$ sample not accepted
8. Save the last accepted sample: $\begin{array}{l}\mathbf{E}_{all}(i,:)=\mathbf{E}_a\end{array}$
Analyze the accepted samples

Examples and Results

To test the Bayesian updating method two examples were considered. The first example is a simple 2DOF system where the basic ideas of the Bayesian updating method could be investigated and well understood. As a second example more complex system was analyzed, a beam discretized using the finite element method.

It was not possible to obtain real measured eigenfrequencies and mode shapes so artificially generated modal data were used instead.

2DOF System

The system is depicted in figure below [3].

The system matrices are given as follows:

$\begin{array}{l}\mathbf{M}=\begin{bmatrix} k_1+k_2 & -k_2 \\ -k_2 & k_2 \end{bmatrix} , \ \mathbf{M}=\begin{bmatrix} m_1 & 0 \\ 0 & m_2 \end{bmatrix} .\end{array}$

For simplicity, only the spring stiffness $\begin{array}{l}k_1\end{array}$ and $\begin{array}{l}k_2\end{array}$ were chosen as random variables and the masses were considered to be deterministic. The prior distribution is considered as in eq. (4) with mean vector set to $\begin{array}{l}\pmb{\mu}=[\mu_{k_1} \ \mu_{k_2}]^T\end{array}$ and coefficients $\begin{array}{l}\pmb{\alpha} = [\alpha_1 \ \alpha_2]^T\end{array}$ . The likelihood function was taken as in eq. (7) with coefficient $\begin{array}{l}\beta\end{array}$ as a parameter.

The artificial measured data were generated using modified inputs $\begin{array}{l}\bar{k}_1\end{array}$ and $\begin{array}{l}\bar{k}_2\end{array}$ for the system matrices, from which the eigenvalues and eigenvectors were computed. The MCMC sampling was performed with number of samples $\begin{array}{l}N_s\end{array}$ . The specific values for the inputs are summarized in the following table and the prior distribution is shown.

$\begin{array}{l}\mu_{k_1}\end{array}$	$\begin{array}{l}100~ [\mathrm{N/m}]\end{array}$
$\begin{array}{l}\mu_{k_2}\end{array}$	$\begin{array}{l}100~ [\mathrm{N/m}]\end{array}$
$\begin{array}{l}\alpha_1\end{array}$	$\begin{array}{l}0.01~ [\mathrm{m}^2/\mathrm{N}^2]\end{array}$
$\begin{array}{l}\alpha_2\end{array}$	$\begin{array}{l}0.01~ [\mathrm{m}^2/\mathrm{N}^2]\end{array}$
$\begin{array}{l}m_1\end{array}$	$\begin{array}{l}100 ~ [\mathrm{kg}]\end{array}$
$\begin{array}{l}m_2\end{array}$	$\begin{array}{l}100 ~ [\mathrm{kg}]\end{array}$
$\begin{array}{l}\beta\end{array}$	$\begin{array}{l}1\end{array}$
$\begin{array}{l}\bar{k}_1\end{array}$	$\begin{array}{l}1.3 \cdot \mu_{k_1} = 130~ [\mathrm{N/m}]\end{array}$
$\begin{array}{l}\bar{k}_2\end{array}$	$\begin{array}{l}1.3 \cdot \mu_{k_2} = 130~ [\mathrm{N/m}]\end{array}$
$\begin{array}{l}N_s\end{array}$	$\begin{array}{l}10^5\end{array}$

In the following figure the histograms of saved parameters $\begin{array}{l}k_1\end{array}$ and $\begin{array}{l}k_2\end{array}$ are plotted.

In the table bellow the posterior mean vector estimate $\begin{array}{l}\pmb{\mu}_p=[\mu_{p,k_1} \ \mu_{p,k_2}]^T\end{array}$ and standard deviation estimates $\begin{array}{l}\sigma_{p,k_1}\end{array}$ and $\begin{array}{l}\sigma_{p,k_2}\end{array}$ are compared with the mean and standard deviation of the prior distribution. Furthermore, the values of the modified input stiffnesses used for the generation of the artificial data are shown.

Prior		Data		Posterior
$\begin{array}{l}\mu_{k_1}\end{array}$	$\begin{array}{l}100~ [\mathrm{N/m}]\end{array}$	$\begin{array}{l}\bar{k}_1\end{array}$	$\begin{array}{l}130~ [\mathrm{N/m}]\end{array}$	$\begin{array}{l}\mu_{p,k_1}\end{array}$	$\begin{array}{l}121.0035~ [\mathrm{N/m}]\end{array}$
$\begin{array}{l}\mu_{k_2}\end{array}$	$\begin{array}{l}100~ [\mathrm{N/m}]\end{array}$	$\begin{array}{l}\bar{k}_2\end{array}$	$\begin{array}{l}130~ [\mathrm{N/m}]\end{array}$	$\begin{array}{l}\mu_{p,k_2}\end{array}$	$\begin{array}{l}128.6601~ [\mathrm{N/m}]\end{array}$
$\begin{array}{l}\sigma_{k_1}\end{array}$	$\begin{array}{l}10~ [\mathrm{N/m}]\end{array}$			$\begin{array}{l}\sigma_{p,k_1}\end{array}$	$\begin{array}{l}6.3482~ [\mathrm{N/m}]\end{array}$
$\begin{array}{l}\sigma_{k_2}\end{array}$	$\begin{array}{l}10~ [\mathrm{N/m}]\end{array}$			$\begin{array}{l}\sigma_{p,k_2}\end{array}$	$\begin{array}{l}3.6250~ [\mathrm{N/m}]\end{array}$

One can notice how the distribution of the parameters was updated based on the measured data. How much does the posterior distribution differ from the prior one strongly depends on the coefficients $\begin{array}{l}\alpha_i\end{array}$ and $\begin{array}{l}\beta\end{array}$ .

Parameter Study

In the following section an influence of both parameters on the posterior mean vector and standard deviation estimates is studied. For convenience, only one parameter $\begin{array}{l}\alpha\end{array}$ was considered for both $\begin{array}{l}k_1\end{array}$ and $\begin{array}{l}k_2\end{array}$ , therefore it holds $\begin{array}{l}\alpha_i=\alpha\end{array}$ .

In figure below, coefficient $\begin{array}{l}\alpha\end{array}$ varies in range $\begin{array}{l}(0.001,5)~ [\mathrm{m}^2/\mathrm{N}^2]\end{array}$ while $\begin{array}{l}\beta=1\end{array}$ is kept constant.

As described before, $\begin{array}{l}\alpha\end{array}$ represents the inverse of the variance of the prior distribution, in other words, it indicates how much trust is put into the prior belief about the individual parameters. Small values of $\begin{array}{l}\alpha\end{array}$ mean large variance of the prior distribution and therefore the parameters are updated more significantly. This is also clear from the figure above, for $\begin{array}{l}\alpha=0.001~ [\mathrm{m}^2/\mathrm{N}^2]\end{array}$ the mean estimates of both stiffnesses are very close to modified inputs $\begin{array}{l}\bar{k}_1\end{array}$ and $\begin{array}{l}\bar{k}_2\end{array}$ based on which the artificial measure modal data were computed. On the other hand, for large values of $\begin{array}{l}\alpha\end{array}$ the prior variance is small therefore only small update is allowed. A similar trend can be observed for the standard deviations.

Now an impact of $\begin{array}{l}\beta\end{array}$ will be examined. It expresses the influence of the likelihood function on the posterior distribution. In the next figure, coefficient $\begin{array}{l}\beta\end{array}$ changes its value in interval $\begin{array}{l}(0.01,10)\end{array}$ and $\begin{array}{l}\alpha=0.1~ [\mathrm{m}^2/\mathrm{N}^2]\end{array}$ .

As can be seen, for increasing $\begin{array}{l}\beta\end{array}$ the posterior mean estimates are approaching the stiffness values corresponding to the measured data. This can also be explained based on the eq. (7), where $\begin{array}{l}\beta\end{array}$ multiplies term $\begin{array}{l}-\[\sum_{i=1}^{k} \sum_{j=1}^{k} \varepsilon_{ij}^2 \]\end{array}$ . This way $\begin{array}{l}\beta\end{array}$ can be seen as a penalization factor. However, for large values of $\begin{array}{l}\beta\end{array}$ numerical instability can occur because the likelihood function values will decrease to zero very rapidly and the computer precision is limited.

Lognormal Prior

As discussed in section Prior Probability Distribution Function, a lognormal prior distribution might be more suitable for description of the parameters $\begin{array}{l}k_1\end{array}$ and $\begin{array}{l}k_2\end{array}$ , since stiffness can only take positive values. In the following, previously investigated system was computed with lognormal prior taken from eq. (5) and eq. (6). Specific values of the distribution parameters are stated in the following table and the related prior distribution is plotted and compared to the normal prior.

$\begin{array}{l}\mu_{k_1}\end{array}$	$\begin{array}{l}100~ [\mathrm{N/m}]\end{array}$	$\begin{array}{l}\mu_{lnk_1}\end{array}$	$\begin{array}{l}4.6002~ [\mathrm{N/m}]\end{array}$
$\begin{array}{l}\mu_{k_2}\end{array}$	$\begin{array}{l}100~ [\mathrm{N/m}]\end{array}$	$\begin{array}{l}\mu_{lnk_2}\end{array}$	$\begin{array}{l}4.6002~ [\mathrm{N/m}]\end{array}$
$\begin{array}{l}\sigma_{k_1}\end{array}$	$\begin{array}{l}10~ [\mathrm{N/m}]\end{array}$	$\begin{array}{l}\sigma_{lnk_1}\end{array}$	$\begin{array}{l}0.0998~ [\mathrm{N/m}]\end{array}$
$\begin{array}{l}\sigma_{k_1}\end{array}$	$\begin{array}{l}10~ [\mathrm{N/m}]\end{array}$	$\begin{array}{l}\sigma_{lnk_1}\end{array}$	$\begin{array}{l}0.0998~ [\mathrm{N/m}]\end{array}$

In figure below histograms of the sampled parameters $\begin{array}{l}k_1\end{array}$ and $\begin{array}{l}k_2\end{array}$ are plotted.

A comparison of the results with different prior distributions is provided in the following table. The deviation of the results is insignificant since the difference between the prior distributions for the specific parameters is minor.

Posterior distribution parameters
	Gaussian prior	Lognormal prior
$\begin{array}{l}\mu_{p,k_1}\end{array}$	$\begin{array}{l}121.0035~ [\mathrm{N/m}]\end{array}$	$\begin{array}{l}121.9600~ [\mathrm{N/m}]\end{array}$
$\begin{array}{l}\mu_{p,k_2}\end{array}$	$\begin{array}{l}128.6601~ [\mathrm{N/m}]\end{array}$	$\begin{array}{l}129.4454~ [\mathrm{N/m}]\end{array}$
$\begin{array}{l}\sigma_{p,k_1}\end{array}$	$\begin{array}{l}6.3482~ [\mathrm{N/m}]\end{array}$	$\begin{array}{l}7.2798~ [\mathrm{N/m}]\end{array}$
$\begin{array}{l}\sigma_{p,k_2}\end{array}$	$\begin{array}{l}3.6250~ [\mathrm{N/m}]\end{array}$	$\begin{array}{l}3.9279~ [\mathrm{N/m}]\end{array}$

Beam example

The Bayesian update using MCMC can also be applied to more complex examples, like a beam with defects. The analyzed case is a beam of length $\begin{array}{l}L\end{array}$ fixed in both ends with a mass density $\begin{array}{l}\rho\end{array}$ , Young's modulus $\begin{array}{l}E_Y\end{array}$ and a rectangular cross-sectional area $\begin{array}{l}A_{acc}\end{array}$ . However, due to a manufacturing defect or holes drill in it, the cross-sectional area $\begin{array}{l}A_{nacc}\end{array}$ of a part of the beam is not accurately known. This case can be solved with a Finite Element (FE) model with a discretization like the one in the image.

In this case, the beam has been discretized by $\begin{array}{l}n=11\end{array}$ elements, of which the white ones have the same area $\begin{array}{l}A_{acc}\end{array}$ and the red ones $\begin{array}{l}A_{nacc}\end{array}$ . The modal properties of interest are the first 4 eigenfrequencies. Without model updating, the "real" modes cannot be accurately known beforehand, so this technique allows a better approach to the problem. All the variables will be modelled with normal distributions of a set mean $\begin{array}{l}\mu\end{array}$ and a standard deviation $\begin{array}{l}\sigma\end{array}$ . The theoretical sides of the cross section are, based in the values used by Marwala in [1], $\begin{array}{l}a=0.0292\end{array}$ m and $\begin{array}{l}b=0.0096\end{array}$ m, and the length of the beam $\begin{array}{l}L=1.1\end{array}$ m. The set of updating parameters is the following:

E	Variable	Mean	Standard deviation	Units
$\begin{array}{l}E_1\end{array}$	$\begin{array}{l}A_{acc}\end{array}$	$\begin{array}{l}2.8032\times10^{-4}\end{array}$	$\begin{array}{l}A_{acc}/20=4.4016\times10^{-5}\end{array}$	$\begin{array}{l}m^2\end{array}$
$\begin{array}{l}E_2\end{array}$	$\begin{array}{l}E_Y\end{array}$	$\begin{array}{l}700\times10^8\end{array}$	$\begin{array}{l}20\times10^8\end{array}$	$\begin{array}{l}N/m^3\end{array}$
$\begin{array}{l}E_3\end{array}$	$\begin{array}{l}\rho\end{array}$	$\begin{array}{l}2700\end{array}$	$\begin{array}{l}25\end{array}$	$\begin{array}{l}kg/m^3\end{array}$
$\begin{array}{l}E_4\end{array}$	$\begin{array}{l}A_{nacc}\end{array}$	$\begin{array}{l}2.8032\times10^{-4}\end{array}$	$\begin{array}{l}A_{nacc}/20=4.4016\times10^{-5}\end{array}$	$\begin{array}{l}m^2\end{array}$

Set-up of the problem

A random sample from this model is taken as the "true" set of parameters. This sample could also have been chosen arbitrarily, as in the previous example. One possible set of variables drawn is, for example:

E	Variable	Value	Units
$\begin{array}{l}E_1\end{array}$	$\begin{array}{l}A_{acc}\end{array}$	$\begin{array}{l}2.8770\times10^{-4}\end{array}$	$\begin{array}{l}m^2\end{array}$
$\begin{array}{l}E_2\end{array}$	$\begin{array}{l}E_Y\end{array}$	$\begin{array}{l}719.15\times10^8\end{array}$	$\begin{array}{l}N/m^3\end{array}$
$\begin{array}{l}E_3\end{array}$	$\begin{array}{l}\rho\end{array}$	$\begin{array}{l}2687.7\end{array}$	$\begin{array}{l}kg/m^3\end{array}$
$\begin{array}{l}E_4\end{array}$	$\begin{array}{l}A_{nacc}\end{array}$	$\begin{array}{l}2.8032\times10^{-4}\end{array}$	$\begin{array}{l}m^2\end{array}$

For this set of parameters, a FEM model of the beam (modelled in Matlab® with the package provided for the course of Modelling and Simulation in Structural Mechanics) gives the following natural frequencies $\begin{array}{l}\omega\end{array}$ and its associated eigenmodes with the plotted shapes:

Eigenfrequencies	[rad/s]
$\begin{array}{l}\omega_1\end{array}$	795
$\begin{array}{l}\omega_2\end{array}$	2191
$\begin{array}{l}\omega_3\end{array}$	4289
$\begin{array}{l}\omega_4\end{array}$	7085

This will be the "truth" against which the full MCMC procedure will be measured. For setting the MCMC, first the $\begin{array}{l}\alpha\end{array}$ and $\begin{array}{l}\beta\end{array}$ parameters have to be chosen. In this case, the chosen parameters are:

Parameter	Variable	Value
$\begin{array}{l}\alpha_1\end{array}$	$\begin{array}{l}A_{acc}\end{array}$	100
$\begin{array}{l}\alpha_2\end{array}$	$\begin{array}{l}\rho\end{array}$	100
$\begin{array}{l}\alpha_3\end{array}$	$\begin{array}{l}E_Y\end{array}$	100
$\begin{array}{l}\alpha_4\end{array}$	$\begin{array}{l}A_{nacc}\end{array}$	30
$\begin{array}{l}\beta\end{array}$	$\begin{array}{l}\varepsilon_{ij}\end{array}$	15

These parameters were chosen in such a way that the probabilities $\begin{array}{l}P(\mathbf{E})\end{array}$ and $\begin{array}{l}P(\mathbf{D}|\mathbf{E})\end{array}$ have the same order of magnitude and the rejection rate is reasonable. Also, $\begin{array}{l}\alpha_4\end{array}$ is chosen significantly smaller than the other $\begin{array}{l}\alpha\end{array}$ to update the defective areas more often. However, this is not sufficient for the convergence of the MCMC to a reasonable value. It was observed, that the MCMC would have a tendency to update the Young's modulus more often than the other variables and to set the cross-sectional areas to 0. This was thought to be related to the big difference between the orders of magnitudes of the parameters. To reduce this effect, a normalization of the error was introduced. The prior probability function is modified to

(13)	$\begin{array}{l}\displaystyle P(\mathbf{E}) = \mathrm{exp} \left[ -\frac{1}{2} \[\sum_{i=1}^{k} \alpha_i \left(\dfrac{E_i - \mu_i}{\mu_i}\right)^2 \] \right].\end{array}$

This is equivalent to dividing every $\begin{array}{l}\alpha_i\end{array}$ by the respective $\begin{array}{l}\mu_i\end{array}$ of the theoretical distribution of the parameters. In case of applying the standard equation, this normalization could have been included in the $\begin{array}{l}\alpha\end{array}$ parameters. For the likelihood probability distribution, we operate analogously:

(14)	$\begin{array}{l}\displaystyle P(\mathbf{D}\|\mathbf{E}) = \mathrm{exp} \left[ - \beta \[\sum_{i=1}^{k} \sum_{j=1}^{k} \left({\dfrac{\left( - \omega_{m,i}^2 \mathbf{M} + \mathbf{K} \right) \pmb{\phi}_{m,i} }{\omega_{m,i}^2}}\right)_j^2 \] \right] ,\end{array}$

In this case, the influence of every mode is normalized with respect to its eigenfrequency. Once more, it is equivalent to apply different $\begin{array}{l}\beta_i\end{array}$ to each of the modes, and varying this factor in each iteration to adjust to the new modes. This ensures that the criteria for updating is uniform for every mode, although it could be modified in case any of the modes would be of special interest. This normalization factor varies with each iteration, but the variance in the eigenfrequencies between two steps is sufficiently enough that the convergence of the MCMC is not compromised. Also, it does not generate a bias in the solution, as the ideal 0 remains in the case when the eigenvalue problem and the parameters are 0, i.e. the modes are the same as the measured ones.

Analysis of results

The MCMC is run for 50000 iterations with a rejection rate of 19%. The results of this simulation are represented in the following plots. In all the line plots, the truth is represented in red, the mean with an interval of confidence of the 95% is represented in blue (continuous for the mean and dashed for the confidence interval) and the random walk in green. First, we analyze the convergence of the parameters vector $\begin{array}{l}\mathbf{E}\end{array}$ :

From the random walk, we can see that the lack of burn out period stops having a significant influence after the first 10000 samples. Also, the areas are modified more often than the other parameters. This is probably due to the greater standard deviation with respect to the mean of their respective distribution. This is in line with a possible real model, as the density and the Young's modulus of the material are usually more accurate than the areas of the beam. If we focus only on the last steps, we get the following plots:

Here, the effect of the parameter $\begin{array}{l}\alpha\end{array}$ is observable. The variation of the least accurate cross-sectional $\begin{array}{l}E_4\end{array}$ is greater than that of $\begin{array}{l}E_1\end{array}$ due to having a smaller $\begin{array}{l}\alpha\end{array}$ . Also, in these plots the convergence of the parameters can be appreciated. If we plot the relative error against our ground truth, we obtain the following plots:

In these plots we see the general convergence of all the parameters. This is especially true for the Young's modulus and the density of the beam, converging against the true measured value. The areas still have room to improve the accuracy, however they are under the 5%, which is acceptable. Taken into account that they don't necessarily need to converge to the original ones, so the model can show enough of a convergence.

If we plot the histograms of the parameters, we can see that the distributions have the expected normal shapes. The measured value is indicated in red once more. Summing up, the results of the sampled set can be gathered in the following table:

E	Truth	Standard deviation (theoretical)	Mean (sampled)	Standard deviation (sampled)	Relative error
$\begin{array}{l}E_1\end{array}$	$\begin{array}{l}2.8770\times10^{-4}\end{array}$	$\begin{array}{l}4.4016\times10^{-5}\end{array}$	$\begin{array}{l}2.8078\times10^{-4}\end{array}$	$\begin{array}{l}3.1318\times10^{-5}\end{array}$	$\begin{array}{l}-2.405\%\end{array}$
$\begin{array}{l}E_2\end{array}$	$\begin{array}{l}719.15\times10^8\end{array}$	$\begin{array}{l}20\times10^8\end{array}$	$\begin{array}{l}717.72\times10^8\end{array}$	$\begin{array}{l}189.63\times10^8\end{array}$	$\begin{array}{l}-0.199\%\end{array}$
$\begin{array}{l}E_3\end{array}$	$\begin{array}{l}2687.7\end{array}$	$\begin{array}{l}25\end{array}$	$\begin{array}{l}2685.6\end{array}$	$\begin{array}{l}317.5\end{array}$	$\begin{array}{l}-0.078\%\end{array}$
$\begin{array}{l}E_4\end{array}$	$\begin{array}{l}2.8032\times10^{-4}\end{array}$	$\begin{array}{l}4.4016\times10^{-5}\end{array}$	$\begin{array}{l}2.7823\times10^{-4}\end{array}$	$\begin{array}{l}5.5\times10^{-5}\end{array}$	$\begin{array}{l}-3.679\%\end{array}$

Surprisingly, the parameters with the biggest precision with respect to the mean value of the distribution have standard deviations greater than the theoretical. This can be directly related, as it can be an indicator that the random walk is trying more extreme values that, together with the smaller variation on of the path, help the convergence to the mean value. With more samples, the standard deviation of the variables will approximate to the theoretical one. An analysis of the correlation of the parameters can be performed.

The scatter plot just shows in this case that the variables seem to be independent, with its dispersion related to the standard deviation of the variables itself.

However, it makes no sense to know the true parameters of the model beforehand and still needing a model updating. Usually, only the modal properties of the system are known, so it is relevant to repeat the same analysis with the first four eigenfrequencies.

In the case of the eigenfrequencies, the convergence seems to be faster, and the random walk appears to be truly random. Also, all the eigenfrequencies are modified roughly in the same way, as expected. The relative error plots show the same trend:

As it can be observed, the relative error is oscillating around the 0%, showing that the updating of the model regarding the eigenfrequency is already converged. This is confirmed when analyzing the contributions of the prior and the likelihood probabilities to the total error. in the samples were the model is already almost converged, the likelihood component of the error vanishes.

The histograms show the same normal distributions. They are situated around the truth one, indicated with the red line. This is confirmed with the comparison with the ground truth:

Eigenfrequency	Truth	Sampled	Standard deviation	Relative error
$\begin{array}{l}\omega_1 [rad/s]\end{array}$	795	804	117	1.212%
$\begin{array}{l}\omega_2[rad/s]\end{array}$	2191	2213	322	1.019%
$\begin{array}{l}\omega_3[rad/s]\end{array}$	4289	4342	626	1.232%
$\begin{array}{l}\omega_4[rad/s]\end{array}$	7085	7163	1034	1.103%

All the values for the eigenfrequencies are within reason. The error is only slightly over 1%, which is a good approximation for the eigenfrequency. The standard deviation is also reasonable, being around a 15% of the mean. Therefore, the updated model can predict accurately the measure eigenfrequencies. However, we must take into account that a finite set of modes does not answer to a unique set of parameters. Different sets can generate the same solution for the eigenvalue problem in a number of modes smaller than the number of degrees of freedom. That explains the variability observed in the posterior distribution of the parameters while the posterior distribution of the eigenfrequencies remains accurate.

The scatter plots reflect the linear dependency of the eigenvalues. However, it is remarkable that the first eigenfrequencies are more independent of each other. Updatings in the parameters seem to affect more the first two eigenfrequencies.

In any case, this example shows the applicability of the Bayesian updating using the Monte Carlo Markov Chain method to real models like a FEM beam. The results can be improved with further tuning of the parameters and a larger number of samples. The results obtained are satisfactory regardless.

References

[1] Marwala, Tshilidzi & Sibisi, Sibusiso. (2005). Finite Element Model Updating Using Bayesian Framework and Modal Properties. Journal of Aircraft. 42. 10.2514/1.11841.

[2] Papaioannou I., Stochastic Finite Element Methods, Lecture notes, Engineering Risk Analysis Group, Technical University of Munich, 2020

[3] Dumond, Patrick & Baddour, Natalie. (2014). A structured approach to design-for-frequency problems using the Cayley-Hamilton theorem. SpringerPlus. 3. 10.1186/2193-1801-3-272.

[4] Jin, Seung-Seop & Ju, Heekun & Jung, Hyung-Jo. (2019). Adaptive Markov chain Monte Carlo algorithms for Bayesian inference: recent advances and comparative study. Structure and Infrastructure Engineering. 10.1080/15732479.2019.1628077.

[5] Roberts, G. O.; Gelman, A.; Gilks, W. R. (1997) Weak convergence and optimal scaling of random walk Metropolis algorithms. Ann. Appl. Probab. 7 , no. 1, 110-120. doi:10.1214/aoap/1034625254. https://projecteuclid.org/euclid.aoap/1034625254

Seitenhierarchie

Bayesian updating using modal data