This function runs the Bayesian Moving Particles algorithm for estimating extreme probability and quantile.
BMP(
dimension,
lsf,
q,
N = 1000,
N.final = N,
N.iter = 30,
adaptive = FALSE,
N.DoE = 5 * dimension,
firstDoE = "uniform",
radius = qnorm(1e-10, lower.tail = FALSE),
X,
y,
covariance = NULL,
learn_each_train = Inf,
km.param = list(nugget.estim = TRUE, multistart = 1, optim.method = "BFGS", coef.trend
= q),
burnin = 20,
fast = TRUE,
sur = list(integrated = TRUE, r = 1, approx.pnorm = FALSE),
lower.tail = TRUE,
save.dir,
plot = FALSE,
plot.lsf = TRUE,
plot.lab = c("x_1", "x_2"),
chi2 = FALSE,
verbose = 1,
breaks
)An object of class list containing the outputs described below:
the estimated conditional expectation of the probability.
the sequence of estimated alpha during the refinement step.
an estimate of the squarred coefficient of variation of alpha.
the sequence of the estimated coefficients of variations.
the sequence of the estimated upper bound of the conditional variance divided by estimated alpha.
the sequence of the estimated integrated h.
a list containing the the sequence of corresponding thresholds and -log probability of the sample minimising the SUR criterion.
a list containing at each iterations number of points tried for the SUR criterion as well as the computational spent.
the reference quantile for the probability estimate.
the empirical cdf, i.e. the estimation of the function q -> E(alpha(q)).
the farthest state reached by the random process. Validity range
for the ecdf is then (-Inf, L_max] or [L_max, Inf).
the last Poisson process generated with N.final particles.
the metamodel approximation of the lsf. A call output is a
list containing the value and the standard deviation.
the final metamodel. An S4 object from DiceKriging. Note
that the algorithm enforces the problem to be the estimation of P[lsf(X)>q]
and so using ‘predict’ with this object will return inverse values if
lower.tail==TRUE; in this scope prefer using directly meta_fun which
handles this possible issue.
the first metamodel with the intial DoE.
a 95% confidence intervalle on the estimate of alpha.
a vector containing the number of moves for each one of the N.batch particles.
the output of the chisq.test function.
the dimension of the input space.
the function defining the RV of interest Y = lsf(X).
a given quantile to estimate the corresponding probability.
the total number of Poisson processes during the refinement step.
the total number of Poisson processes for the final alpha estimate.
the total number of iteration of the algorithm, ie that total number of
calls to the lsf will be N.DoE + N.iter*r.
if the algorithm should stop automatically if the stopping criterion is verified, precisely the mean probability of misclassification of the particles being over a given threshold.
the number of points for the initial Design of Experiment
default is "uniform" for a random
uniform sampling over a sphere of radius radius. Also available "maximim" for a maximim LHS.
the size of the radius of the sphere for uniform DoE or the semi length of the interval on each dimension for maximin LHS
(optional) a first Design of Experiemnt to be used instead of building a new DoE
the value of lsf on the X
(optional) to give a covariance kernel for the km object.
a integer: after this limit the covariance parameters are not learnt any more and model is just updated with the new datapoints.
(optional) list of parameters to be passed to DiceKriging::km.
a burnin parameter for Markov Chain drawing of the metamodel based
Poisson process (this does not change the number of calls to lsf).
in current implementation it appears that the call to the metamodel is
faster when doing batch computation. This parameter lets do the Markov chain the other way
around: instead of first selecting a starting point and then applying burnin times the
transition kernel, it creates a working population by apply the kernel to all the
particles and then makes some moves with the generated discretised distribution.
a list containing any parameters to be passed to estimateSUR. Default is
sur$integrated=TRUE and sur$r=1 for a one step ahead integrated SUR criterion.
as for pxxxx functions, TRUE for estimating P(lsf(X) < q), FALSE for P(lsf(X) > q).
(optional) a directory to save the X and y at each iteration.
to plot the DoE and the updated model.
to plot the contour of the true lsf. Note that this requires its
evaluation on a grid and should be used only on toy examples.
the labels of the axis for the plot.
for a chi2 test on the number of events.
controls the level of outputs of the algorithm.
optional, for the final histogram if chi2 == TRUE.
Clement WALTER clementwalter@icloud.com
The Bayesian Moving Particles algorithm uses the point process framework for rare event to iteratively estimate the conditional expectation of the (random) limit-state function, to quantify the quality of the learning and to propose a new point to be added to the model with a SUR criterion.
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.
C. Walter:
Moving Particles: a parallel optimal Multilevel Splitting
method with application in quantiles estimation and meta-model
based algorithms
Structural Safety, 55, 10-25.
J. Bect, L. Li and E. Vazquez:
Bayesian subset simulation
arXiv preprint arXiv:1601.02557
SubsetSimulation
MonteCarlo
IRW
MP
# Estimate P(g(X)<0)
if (FALSE) p <- BMP(dimension = 2, lsf = kiureghian, q = 0, N = 100, N.iter = 30, plot = TRUE)
# More extreme event
if (FALSE) p <- BMP(dimension = 2, lsf = waarts, q = -4, N = 100, N.iter = 50, plot = TRUE)
# One can also estimate probability of the form P(g(X)>q)
if (FALSE) p <- BMP(dimension = 2, lsf = cantilever, q = 1/325, N = 100, N.iter = 30, plot = TRUE)
Run the code above in your browser using DataLab