# 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 () {
logn(x1,mode=1.3,sd=1.),
logn(x2,mode=0.8,sd=1.)
};

### 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));
}::{k1,k2,b,c,f1,f2,J});

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

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
};

statsmp(is_x1);

#### CDF of x1

Fesslix parameter file
echo "
********************************
* Plot CDF of x1               *
********************************
";

# First, we need to sort the samples in 'is_x1':
sortsmp(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
smpplot(
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