mistral (version 1.1-1)

SubsetSimulation: Subset Simulation Monte-Carlo

Description

Estimate a failure probability with a Monte-Carlo method applied on nested subsets.

Usage

SubsetSimulation(dimension,
                 limit_state_function,
                 proposal_pdf,
                 pdf            = dnorm,
                 rpdf           = rnorm,
                 cutoff_prob    = 0.1,
                 n_init_samples = 10000,
				         burnin 	     	= 20,
				         thinning 		  = 4,
                 plot           = FALSE,
                 output_dir     = NULL,
			        	 verbose 		    = 0)

Arguments

dimension

an integer giving the dimension of the input space.

limit_state_function

the failure function.

proposal_pdf

proposal PDF for Metropolis-Hastings algorithm. Default random walk is uniform on a radius of 2.

pdf

PDF of the input space to be used in Metropolis-Hastings algorithm.

rpdf

a random generator for the input space PDF.

cutoff_prob

the cut-off probability for each subset.

n_init_samples

number of samples to be used for each subset.

burnin

burnin parameter for Metropolis-Hastings algorithm.

thinning

thinning parameter for Metropolis-Hastings algorithm.

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

output_dir

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

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:
• proba
  The estimated failure probability.
• cov
  The coefficient of variation of the Monte-Carlo probability estimate.
• Ncall
  The total number of calls to the limit_state_function.

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 cutoff_prob. 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)

See Also

MonteCarlo S2MART

Examples

#Try Subset Simulation Monte Carlo on a given function and change number of points.
#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 = list()
res[[1]] = SubsetSimulation(2,kiureghian,n_init_samples=10000)
res[[2]] = SubsetSimulation(2,kiureghian,n_init_samples=100000)
res[[3]] = SubsetSimulation(2,kiureghian,n_init_samples=500000)

#Try Subset Simulation Monte Carlo on a given function with different failure level
#Limit state function defined by Waarts :
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[[1]] = SubsetSimulation(2,waarts,failure=0,plot=TRUE)
res[[2]] = SubsetSimulation(2,waarts,failure=1,plot=TRUE)
res[[3]] = SubsetSimulation(2,waarts,failure=-1,plot=TRUE)