MonteCarlo: Crude Monte Carlo method

In clemlaflemme/mistral: Methods in Structural Reliability

Description

Estimate a failure probability using a crude Monte Carlo method.

Usage

MonteCarlo(
  dimension,
  lsf,
  N_max = 5e+05,
  N_batch = foreach::getDoParWorkers(),
  q = 0,
  lower.tail = TRUE,
  precision = 0.05,
  plot = FALSE,
  output_dir = NULL,
  save.X = TRUE,
  verbose = 0
)

Arguments

dimension	the dimension of the input space.
lsf	the function defining safety/failure domain.
N_max	maximum number of calls to the lsf.
N_batch	number of points evaluated at each iteration.
q	the quantile.
lower.tail	as for pxxxx functions, TRUE for estimating P(lsf(X) < q), FALSE for P(lsf(X) > q).
precision	a targeted maximum value for the coefficient of variation.
plot	to plot the contour of the lsf as well as the generated samples.
output_dir	to save a copy of the plot in a pdf. This name will be pasted with "_Monte_Carlo_brut.pdf".
save.X	to save all the samples generated as a matrix. Can be set to FALSE to reduce output size.
verbose	to control the level of outputs in the console; either 0 or 1 or 2 for almost no outputs to a high level output.

Details

This implementation of the crude Monte Carlo method works with evaluating batchs of points sequentialy until a given precision is reached on the final estimator

Value

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

p	the estimated probabilty.
ecdf	the empiracal cdf got with the generated samples.
cov	the coefficient of variation of the Monte Carlo estimator.
Ncall	the total numnber of calls to the lsf, ie the total number of generated samples.
X	the generated samples.
Y	the value lsf(X).

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.

Author(s)

Clement WALTER [e-mail]

References

• R. Rubinstein and D. Kroese:
  Simulation and the Monte Carlo method
  Wiley (2008)
  

See Also

SubsetSimulation foreach

Examples

#First some considerations on the usage of the lsf. 
#Limit state function defined by Kiureghian & Dakessian :
# Remember you have to consider the fact that the input will be a matrix ncol >= 1
lsf_wrong = function(x, b=5, kappa=0.5, e=0.1) {
  b - x[2] - kappa*(x[1]-e)^2 # work only with a vector of lenght 2
}
lsf_correct = function(x){
  apply(x, 2, lsf_wrong)
}
lsf = function(x, b=5, kappa=0.5, e=0.1) {
  x = as.matrix(x)
  b - x[2,] - kappa*(x[1,]-e)^2 # vectorial computation, run fast
}

y = lsf(X <- matrix(rnorm(20), 2, 10))
#Compare running time
## Not run: 
  require(microbenchmark)
  X = matrix(rnorm(2e5), 2)
  microbenchmark(lsf(X), lsf_correct(X))

## End(Not run)

#Example of parallel computation
require(doParallel)
lsf_par = function(x){
 foreach(x=iter(X, by='col'), .combine = 'c') %dopar% lsf(x)
}

#Try Naive Monte Carlo on a given function with different failure level
## Not run: 
  res = list()
  res[[1]] = MonteCarlo(2,lsf,q = 0,plot=TRUE)
  res[[2]] = MonteCarlo(2,lsf,q = 1,plot=TRUE)
  res[[3]] = MonteCarlo(2,lsf,q = -1,plot=TRUE)
  

## End(Not run)


#Try Naive Monte Carlo on a given function and change number of points.
## Not run: 
  res = list()
  res[[1]] = MonteCarlo(2,lsf,N_max = 10000)
  res[[2]] = MonteCarlo(2,lsf,N_max = 100000)
  res[[3]] = MonteCarlo(2,lsf,N_max = 500000)

## End(Not run)

clemlaflemme/mistral documentation built on Jan. 3, 2020, 9:13 a.m.