# Subset Simulation Monte Carlo Estimate a probability of failure with the Subset Simulation algorithm (also known as Multilevel Splitting or Sequential Monte Carlo for rare events). ```r 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 ) ``` ## 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 ## Returns 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 ## 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. ## Note Problem is supposed to be defined in the standard space. If not, use `UtoX` to do so. Furthermore, each time a set of vector is defined as a matrix, nrow = `dimension` and ncol = number of vector to be consistent with `as.matrix` transformation of a vector. Algorithm calls lsf(X) (where X is a matrix as defined previously) and expects a vector in return. This allows the user to optimise the computation of a batch of points, either by vectorial computation, or by the use of external codes (optimised C or C++ codes for example) and/or parallel computation; see examples in MonteCarlo . ## Examples ```r #Try Subset Simulation Monte Carlo on a given function and change number of points. ## Not run: res = list() res[[1]] = SubsetSimulation(2,kiureghian,N=10000) res[[2]] = SubsetSimulation(2,kiureghian,N=100000) res[[3]] = SubsetSimulation(2,kiureghian,N=500000) ## End(Not run) # Compare SubsetSimulation with MP ## Not run: 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 ## End(Not run) ``` ## 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` ## Author(s) Clement WALTER clementwalter@icloud.com

SubsetSimulation function