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.
Table of Contents
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
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.)
};
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).
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});
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
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
};
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
};
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);
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
********************************
* 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
