MonteCarlo {mistral}

Crude Monte Carlo method

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)