SubsetSimulation: Subset Simulation Monte Carlo

Description

Estimate a probability of failure with the Subset Simulation algorithm (also known as Multilevel Splitting or Sequential Monte Carlo for rare events).

Usage

SubsetSimulation(
  dimension,
  lsf,
  p_0 = 0.1,
  N = 10000,
  q = 0,
  lower.tail = TRUE,
  K,
  thinning = 20,
  save.all = FALSE,
  plot = FALSE,
  plot.level = 5,
  plot.lsf = TRUE,
  output_dir = NULL,
  plot.lab = c("x", "y"),
  verbose = 0
)

Value

An object of class list containing the failure probability and some more outputs as described below:

p: the estimated failure probability.
cv: the estimated coefficient of variation of the estimate.
Ncall: the total number of calls to the lsf.
X: the working population.
Y: the value lsf(X).
Xtot: if save.list==TRUE, all the Ncall samples generated by the algorithm.
Ytot: the value lsf(Xtot).
sigma.hist: if default kernel is used, sigma is initialized with 0.3 and then further adaptively updated to have an average acceptance rate of 0.3

Arguments

dimension: the dimension of the input space.
lsf: the function defining failure/safety domain.
p_0: a cutoff probability for defining the subsets.
N: the number of samples per subset, ie the population size for the Monte Carlo estimation of each conditional probability.
q: the quantile defining the failure domain.
lower.tail: as for pxxxx functions, TRUE for estimating P(lsf(X) < q), FALSE for P(lsf(X) > q)
K: a transition Kernel for Markov chain drawing in the regeneration step. K(X) should propose a matrix of candidate sample (same dimension as X) on which lsf will be then evaluated and transition accepted of rejected. Default kernel is the one defined K(X) = (X + sigma*W)/sqrt(1 + sigma^2) with W ~ N(0, 1).
thinning: a thinning parameter for the the regeneration step.
save.all: if TRUE, all the samples generated during the algorithms are saved and return at the end. Otherwise only the working population is kept at each iteration.
plot: to plot the generated samples.
plot.level: maximum number of expected levels for color consistency. If number of levels exceeds this value, the color scale will change according to ggplot2 default policy.
plot.lsf: a boolean indicating if the lsf should be added to the plot. This requires the evaluation of the lsf over a grid and consequently should be used only for illustation purposes.
output_dir: to save the plot into a pdf file. This variable will be paster with "_Subset_Simulation.pdf"
plot.lab: the x and y labels for the plot
verbose: Either 0 for almost no output, 1 for medium size output and 2 for all outputs

Author

Clement WALTER clementwalter@icloud.com

Details

This algorithm uses the property of conditional probabilities on nested subsets to calculate a given probability defined by a limit state function.

It operates iteratively on ‘populations’ to estimate the quantile corresponding to a probability of p_0. Then, it generates samples conditionnaly to this threshold, until found threshold be lower than 0.

Finally, the estimate is the product of the conditional probabilities.

References

S.-K. Au, J. L. Beck:
Estimation of small failure probabilities in high dimensions by Subset Simulation
Probabilistic Engineering Mechanics (2001)
A. Guyader, N. Hengartner and E. Matzner-Lober:
Simulation and estimation of extreme quantiles and extreme probabilities
Applied Mathematics and Optimization, 64(2), 171-196.
F. Cerou, P. Del Moral, T. Furon and A. Guyader:
Sequential Monte Carlo for rare event estimation
Statistics and Computing, 22(3), 795-808.

Examples

Run this code

#Try Subset Simulation Monte Carlo on a given function and change number of points.
 
if (FALSE) {
 res = list()
 res[[1]] = SubsetSimulation(2,kiureghian,N=10000)
 res[[2]] = SubsetSimulation(2,kiureghian,N=100000)
 res[[3]] = SubsetSimulation(2,kiureghian,N=500000)
}

# Compare SubsetSimulation with MP
if (FALSE) {
p <- res[[3]]$p # get a reference value for p
p_0 <- 0.1 # the default value recommended by Au and Beck
N_mp <- 100
# to get approxumately the same number of calls to the lsf
N_ss <- ceiling(N_mp*log(p)/log(p_0))
comp <- replicate(50, {
ss <- SubsetSimulation(2, kiureghian, N = N_ss)
mp <- MP(2, kiureghian, N = N_mp, q = 0)
comp <- c(ss$p, mp$p, ss$Ncall, mp$Ncall)
names(comp) = rep(c("SS", "MP"), 2)
comp
})
boxplot(t(comp[1:2,])) # check accuracy
sd.comp <- apply(comp,1,sd)
print(sd.comp[1]/sd.comp[2]) # variance increase in SubsetSimulation compared to MP

colMeans(t(comp[3:4,])) # check similar number of calls
}

Run the code above in your browser using DataLab

Last chance! 50% off unlimited learning