MetaIS: Metamodel based Impotance Sampling

Description

Estimate failure probability by MetaIS method.

Usage

MetaIS(dimension,
      limit_state_function,
      N                 = 500000,
      N_alpha           = 100,
      N_DOE             = 2*dimension,
      N1                = N_DOE*30,  
      Ru                = 8,
      Nmin              = 30,
      Nmax              = 200,
      Ncall_max         = 1000,
      precision         = 0.05,
      N_seeds           = 1,
      Niter_seed        = 10000,
      N_alphaLOO        = 5000,
      K_alphaLOO        = 2*dimension,
      alpha_int         = c(0.1,10),
      k_margin          = 1.96,
      learn_db          = NULL,
      lsf_value         = NULL,
      failure           = 0,
      meta_model        = NULL,
      kernel            = "matern5_2",
      learn_each_train  = FALSE,
      limit_fun_MH      = NULL,
      sampling_strategy = "MH",
      seeds             = NULL,
      seeds_eval        = NULL,
      burnin            = 20,
      thinning          = 4,
      plot              = FALSE,
      limited_plot      = FALSE,
      add               = FALSE,
      output_dir        = NULL,
      z_MH              = NULL,
      z_lsf             = NULL,
      verbose		        = 0)

Arguments

dimension

an integer giving the dimension of the input space.

limit_state_function

the failure fonction.

an integer defining the size of Monte-Carlo population for augmented failure probability estimation.

N_alpha

an integer defining the initial size of Monte-Carlo population for alpha estimate.

N_DOE

an integer defining the size of initial DOE got by clustering of the N1 samples.

an integer defining the size of the initial uniform population sampled in a hypersphere of radius Ru.

a real defining the radius of the hypersphere for the initial uniform sampling.

Nmin

an integer defining the minimum number of calls to the limit_state_function for the construction step.

Nmax

an integer defining the maximum number of calls to the limit_state_function for the construction step.

Ncall_max

an integer defining the maximum number of calls to the limit_state_function for the whole algorithm.

precision

a real defining the targeted maximal value of the coefficient of variation.

N_seeds

an integer defining the number of seeds for Metropolis-Hastings algorithm while generating into the margin (according to MarginProbability*gauss density).

Niter_seed

maximum number of iterations for the research of a seed for alphaLOO refinement sampling.

N_alphaLOO

an integer defining the number of points to sample at each refinement step.

K_alphaLOO

an integer defining the number of clusters at each refinement step.

alpha_int

a 2-d real vector defining the range for alphaLOO criterion : alpha_int[1] < alphaLOO < alpha_int[2].

k_margin

a real defining the margin width according to standard gaussian law; default is 1.96 which means points outside of the margin are classified with more than 97,5%.

learn_db

optional. A matrix of already known points, with dim : dimension x number_of_vector.

lsf_value

values of the limit state function on the vectors given in learn_db.

failure

the value defining the failure domain F = { x | limit_state_function(x) < failure }.

meta_model

optional. If a kriging based metamodel has already been fitted to the data (from DiceKriging package) it can be given as an input to keep the same parameters.

kernel

a specified kernel to be used for the metamodel. See DiceKriging for available options.

learn_each_train

specify if hyperparameters of the model should be evaluated each time points are added to the learning database ("TRUE") or only the first time ("FALSE").

limit_fun_MH

optional. If the working space is to be reduced to some subset defining by a function, eg. in case of use in a Subset Simulation algorithm. As for the limit_state_function, failure domain is defined by points whom values of limit_fun_MH

sampling_strategy

either "AR" or "MH", to specify which sampling strategy is to be used when generating Monte-Carlo population in a case of subset simulation : "AR" stands for accept-reject while "MH" stands for Metropolis-Hastings.

seeds

optional. If sampling_strategy=="MH", seeds from which MH algorithm starts. This should be a matrix with nrow = dimension and ncol = number of vector.

seeds_eval

optional. The value of the limit_fun_MH on the seeds.

burnin

a burnin parameter for Metropolis-Hastings algorithm. To be used only for Monte-Carlo population sampling and set to 0 elsewhere.

thinning

a thinning parameter for Metropolis-Hastings algorithm. thinning = 0 means no thinning. To be used only for Monte-Carlo population sampling and set to 0 elsewhere.

plot

a boolean parameter specifying if function and samples should be plotted. The plot is refreshed at each iteration with the new data. Note that this option is only to be used when working on light limit state functions as it requires the c

limited_plot

only a final plot with limit_state_function, final DOE and metamodel. Should be used with plot==FALSE.

add

optional. "TRUE" if plots are to be added to the current active device.

output_dir

optional. If plots are to be saved in .jpeg in a given directory. This variable will be pasted with "_MetaIS.jpeg" to get the full output directory.

z_MH

optional. For plots, if metamodel has already been evaluated on the grid then z_MH (from outer function) can be provided to avoid extra computational time.

z_lsf

optional. For plots, if LSF has already been evaluated on the grid then z_lsf (from outer function) can be provided to avoid extra computational time.

verbose

Eiher 0 for an almost no output message, or 1 for medium size or 2 for full size

Value

An object of class list containing the failure probability and some more outputs as described below:
probaThe estimated failure probability.
covThe coefficient of variation of the Monte-Carlo probability estimate.
NcallThe total number of calls to the limit_state_function.
learn_dbThe final learning database, ie. all points where limit_state_function has been calculated.
lsf_valueThe value of the limit_state_function on the learning database.
meta_funThe metamodel approximation of the limit_state_function. A call output is a list containing the value and the standard deviation.
meta_modelThe final metamodel. An S4 object from DiceKriging.
pointsPoints in the failure domain according to the metamodel.
meta_evalEvaluation of the metamodel on these points.
z_metaIf plot==TRUE, the evaluation of the metamodel on the plot grid.

Details

MetaIS is an Important Sampling based probability estimator. It makes use of a kriging surogate to approximate the optimal density function, replacing the indicatrice by its kriging pendant, the probability of being in the failure domain. In this context, the normallizing constant of this quasi-optimal PDF is called the augmented failure probability and the modified probability alpha. After a first uniform Design of Experiments, MetaIS uses an alpha-Leave-One-Out criterion combined with a margin sampling strategy to refine a kriging-based metamodel. Samples are generated according to the weighted margin probability with Metropolis-Hastings algorithm and some are selected by clustering; the N_seeds are got from an accept-reject strategy on a standard population. Once criterion is reached or maximum number of call done, the augmented failure probability is estimated with a crude Monte-Carlo. Then, a new population is generated according to the quasi-optimal instrumenal PDF; burnin and thinning are used here and alpha is evaluated. While the coefficient of variation of alpha estimate is greater than a given threshold and some computation spots still available (defined by Ncall_max) the estimate is refined with extra calculus. The final probability is the product of p_epsilon and alpha, and final squared coefficient of variation is the sum of p_epsilon and alpha one's.

References

V. Dubourg: Meta-modeles adaptatifs pour l'analyse de fiabilite et l'optimisation sous containte fiabiliste PhD Thesis, Universite Blaise Pascal - Clermont II,2011
V. Dubourg, B. Sudret, F. Deheeger: Metamodel-based importance sampling for structural reliability analysis Original Research Article Probabilistic Engineering Mechanics, Volume 33, July 2013, Pages 47-57
V. Dubourg, B. Sudret: Metamodel-based importance sampling for reliability sensitivity analysis. Accepted for publication in Structural Safety, special issue in the honor of Prof. Wilson Tang.(2013)
V. Dubourg, B. Sudret and J.-M. Bourinet: Reliability-based design optimization using kriging surrogates and subset simulation. Struct. Multidisc. Optim.(2011)

Examples

Run this code

#Limit state function defined by Kiureghian & Dakessian :
kiureghian = function(x, b=5, kappa=0.5, e=0.1) {
x = as.matrix(x)
b - x[2,] - kappa*(x[1,]-e)^2
}

res = MetaIS(dimension=2,limit_state_function=kiureghian,plot=TRUE)

#Compare with crude Monte-Carlo reference value
N = 500000
U = matrix(rnorm(dimension*N),dimension,N)
G = apply(U,2,kiureghian)
P = mean(G<0)
cov = sqrt((1-P)/(N*P))

#See impact of kernel choice with Waarts function :
waarts = function(x) {
x = as.matrix(x)
apply(x, 2, function(u) {min(
  	(3+(u[1]-u[2])^2/10 - (u[1]+u[2])/sqrt(2)),
		(3+(u[1]-u[2])^2/10 + (u[1]+u[2])/sqrt(2)),
		u[1]-u[2]+7/sqrt(2),
		u[2]-u[1]+7/sqrt(2))})
}

res = list()
res$matern5_2 = MetaIS(2,waarts,plot=TRUE)
res$matern3_2 = MetaIS(2,waarts,kernel="matern3_2",plot=TRUE)
res$gaussian = MetaIS(2,waarts,kernel="gauss",plot=TRUE)
res$exp = MetaIS(2,waarts,kernel="exp",plot=TRUE)

#Compare with crude Monte-Carlo reference value
N = 500000
U = matrix(rnorm(dimension*N),dimension,N)
G = apply(U,2,waarts)
P = mean(G<0)
cov = sqrt((1-P)/(N*P))

Run the code above in your browser using DataLab