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 [1]. 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

[1] 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)