mistral (version 2.2.2)

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

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.
  

See Also

IRW MP MonteCarlo

Examples

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