Tutorial 2a – Likelihood function formulated directly in Fesslix

In this tutorial we formulate the likelihood function of the problem introduced in Tutorial 2 directly in Fesslix. BUS-SuS is employed to estimate the probability of failure.

1   Step by Step Instruction

1.1   Constants

We start by defining scalar variables:

Fesslix parameter file
const Nsmp = 1e3;  # number of samples used per level in Subset Simulation
const pthr = 10%;  # probability of the intermediate Subset levels

1.2   Prior distribution

The prior distribution of this problem consists of two independent log-Normal random variables, denoted x1 and x2. We group the prior random variables in a set called RBS.

Fesslix parameter file
rbrv_set RBS () {

1.3   Likelihood function

In this tutorial, we define the limit-state function directly within Fesslix.

In the following we use the keyword fun to define a scalar function labeled 'likeli' that takes two parameters (the current realizations of the two-dimensional prior distribution).

Fesslix parameter file
const h1 = 3.13;     # the 1st observed eigenfrequency
const h2 = 9.83;     # the 2nd observed eigenfrequency
const sgma_e = 1/16; # the standard deviation associated with the observation errors
const k = 29.7e6;    # multiplying this value with the random variables x1 and x2
                     # gives the stiffness values
const m1 = 16.531e3; # the lumped mass associated with the 1st DOF
const m2 = 16.131e3; # the lumped mass associated with the 2nd DOF
const a = m1*m2;     # a constant that is needed in the model

fun likeli(2) = objexec(:res,{
    # Evaluate the eigenfrequencies f1 and f2
    # based on the parameters of the function of $1 and $2
      const k1 = $1*k;
      const k2 = $2*k;
      const b = -m1*k2-m2*(k1+k2);
      const c = sqrt(b^2-4*a*k1*k2);
      const f1 = sqrt((-b-c)/(2*a))/2/PI;
      const f2 = sqrt((-b+c)/(2*a))/2/PI;
    # evaluate the Likelihood function
      const J = ((f1/h1)^2-1)^2+((f2/h2)^2-1)^2;
      const res = exp(-J/(2*sgma_e^2));

In the code above, objexec is a function that executes the code in curly brakes before it returns the value of the scalar variable res. The scalar variables k1, k2, b, c, f1, f2, J and res are treated as local variables inside of objexec: Their value is not modified outside of objexec.

1.4   Bayesian inference

We tackle the reliability problem with BUS-SuS. All Bayesian inference methods are contained in the module BayUp which needs to be loaded explicitly:

Fesslix parameter file
loadlib "BayUp"; # load the 'Bayesian Updating'-module

We need to create a set that will contain the posterior distribution after the inference completed. The set is labeled efup and uses the set RBS as prior distribution:

Fesslix parameter file
bayup_new efup {
   rbrvsets = "RBS"; # The set of random variables to update

In a next step, we link the likelihood function to the updating problem.

Fesslix parameter file
bayup_likelihood efup = likeli( rbrv(RBS::x1), rbrv(RBS::x2) );

Finally, we can start the Bayesian inference:

Fesslix parameter file
bayup_update efup: uBUS (
  Nsmp*pthr;            # The number of chains used at each level in Subset Simulation.
  1/pthr;               # The length of each chain.
  Nsmp                  # number of samples in the final conditioning stage
) {
  mcmc_algo = "csusmh"; # use conditional sampling in U-space as MCMC algorithm
  adaptive_meth = log;  # use the 'log' scheme to update the scaling
                        # of the MCMC proposal distribution
  randomize = init;     # the initial ordering of the posterior sampling is randomized
  ext_out = true;       # use extended output

1.5   Work with posterior samples

Compute the mean and standard deviation of the posterior samples of x1

Fesslix parameter file
loadlib "RND"; # in the followin we need functionality from this module

echo "
*Compute mean and std.dev of x1*

ivstream is_x1() {      # create a vector-input stream with name 'is_x1'
  nreserve = Nsmp;      # pre-allocate space for Nsmp entries

sfor(i;Nsmp) {
  rnd_smp "efup";       # retrieve a sample from the posterior set
  ivstream_append is_x1 += rbrv(RBS::x1); # append the sample to the vector-input stream


CDF of x1

Fesslix parameter file
echo "
* Plot CDF of x1               *

# First, we need to sort the samples in 'is_x1':

# Next, we specify where to store the data-file:
  filestream fs1 ( "tutorial_2_cdf_x1.out" );

# Finally, we create the CDF data-file
    is_x1,      # plot this data-set
    type=1      # plot the CDF of this data-set
  ) {
    stream=fs1; # send the output to stream 'fs1'.

# You can plot the results with GNUplot:
# p [][0:1] 'tutorial_2_cdf_x1.out' u 1:2 w l

2   The complete parameter file

fesslix.org – Home  |  Contact  |  Impressum  |  © 2015-2017 Wolfgang Betz