# Multilevel uncertainty quantification with Bayesian inversion

Hello all,

I am trying to solve the Bayesian multilevel model updating approach, which was presented by Nagel and Sudret in . I am applying the example of the paper (Section 8) using the marginalised formulation proposed. The example consists in a simply supported beam, in which variability is in the elastic moduli Ei across the beam (i=1,…,n – where n is the number of experiments). The deflection of the beam is calculated as followed:

Displacement = yi = @(x) F.*x/2./(48*E*I).*(3*L^2-4*(x).^2)


where x is the position along the main axis of the beam (I am considering 3 positions along the beam: 0.25 m, 0.5 m, and 0.75 m). The data of the beam are:

L = 1; h = 0.1; b = 0.1; I = 1/12*b*h^3; F = 30e3;


The Young’s moduli (E_i) are independently sampled from a lognormal distribution with:
\mu_E = 15e9 Pa and \sigma_E = 3e9 Pa.

The experimental data are perturbed with a noise term, which is sampled from a Gaussian distribution with \mu=0 and \sigma = 0.1 mm. The quantity of interest is the hyperparameter \theta_E = (\mu_E,\sigma_E) ((E_i are defined as nuisance variables in the paper).

The first step in applying the multilevel approach is the definition of a suitable prior distribution \pi(\theta_E), in this case, the marginals are independent, so the Bayesian prior is \pi(\theta_E) = \pi(\mu_E)\pi(\sigma_E), these marginals are given in uniform distributions:

\pi(\mu_E)=U[0,100]GPa and \pi(\sigma_E)=U[0,30]GPa

The distributions f(E_i|\theta_E)=LN(E_i|\mu_E,\sigma_E) and \pi(\theta_E) represent the available structural and parametric prior knowledge, respectively. The posterior has the form: \pi(\theta_E|y_i) where y_i are the experimental data (displacement measurements). It can be directly sampled from the marginals of the posterior \pi(<E_i>,\theta_E|<y_i>). The next step is the definition of the posterior density (like for any Bayesian updating), which is marginalised; moreover, thanks to the marginal prior distribution, the unscaled version of the marginal posterior may be written as:

\pi(\theta_E|y_i) = L(\theta_E;<y_i>)\pi(\theta_E)

The next step is to calculate the marginalised likelihood to sample the marginal posterior. Since the marginalised likelihood is a product of integrals, it is not possible to perform the marginalisation analytically, so the stochastic integration via Monte Carlo method is applied.

I am trying to solve the example with the marginalised formulation instead of the joint one because I would like to apply the marginalised formulation of the problem in the application I am working on. However, the standard deviation of the hyperparameter \theta_E is not correctly identified (the solution jumps close to zero). I am not sure why this happens. Instead, the mean is perfectly obtained. Following the example of the paper, I applied Markov Chain Monte Carlo sampling (Metropolis-Hastings algorithm), and the marginalised likelihood is solved with Monte Carlo integration. These are the option for the solution:

Solver.Type = 'MCMC';
Solver.MCMC.Sampler = 'MH';
Solver.MCMC.Steps = 5e3;
Solver.MCMC.NChains = 5;
Solver.MCMC.Proposal.PriorScale = 1e-3;


Here is the solution I obtained,

%------------------- Posterior Marginals
---------------------------------------------------------------
| Parameter | Mean    | Std     | (0.025-0.97) Quant. | Type  |
---------------------------------------------------------------
| meanE     | 1.5e+10 | 6.2e+08 | (1.4e+10 - 1.6e+10) | Model |
| stdE      | 5e+08   | 3.7e+08 | (2.6e+07 - 1.4e+09) | Model |
---------------------------------------------------------------

%------------------- Point estimate
----------------------------------------
| Parameter | Mean    | Parameter Type |
----------------------------------------
| meanE     | 1.5e+10 | Model          |
| stdE      | 5e+08   | Model          |
----------------------------------------

%------------------- Correlation matrix (model parameters)
---------------------------
|       |  meanE    stdE  |
---------------------------
| meanE |  1        0.13  |
| stdE  |  0.13     1     |
---------------------------

%------------------- Convergence
Gelman-Rubin MPSRF:                            1.005291e+00


As I said above, the mean of the hyperparameters is perfectly detected, instead, the standard deviation is much lower than expected. Do you know what I am doing wrong?

I am sorry for this very long post; I hope everything is well explained (I attached the codes (the loglikelihood and the main code) – there is also the solution I obtained since the simulation took some time to run).

Let me know if you need more explanation,

Many thanks in advance,

Irma

 Nagel, Joseph B., and Bruno Sudret. “A unified framework for multilevel uncertainty quantification in Bayesian inverse problems.” Probabilistic Engineering Mechanics 43 (2016): 68-84. Redirecting

Sudret_uq_prova.zip (775.4 KB)