Learn R Programming

mistral (version 2.2.2)

MonteCarlo: 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
)

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

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.

Author

Clement WALTER clementwalter@icloud.com

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

References

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

See Also

SubsetSimulation foreach

Examples

Run this code
#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
if (FALSE) {
  require(microbenchmark)
  X = matrix(rnorm(2e5), 2)
  microbenchmark(lsf(X), lsf_correct(X))
}

#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
if (FALSE) {
  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)
  
}


#Try Naive Monte Carlo on a given function and change number of points.
if (FALSE) {
  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)
}

Run the code above in your browser using DataLab