This function runs the Moving Particles algorithm for estimating extreme probability and quantile.

1 2 3 |

`dimension` |
the dimension of the input space. |

`lsf` |
the function defining the RV of interest Y = lsf(X). |

`N` |
the total number of particles, |

`N.batch` |
the number of parallel batches for the algorithm. Each batch will then
have |

`p` |
a given probability to estimate the corresponding quantile (as in qxxxx functions). |

`q` |
a given quantile to estimate the corresponding probability (as in pxxxx functions). |

`lower.tail` |
as for pxxxx functions, TRUE for estimating P(lsf(X) < q), FALSE for P(lsf(X) > q). |

`Niter_1fold` |
a function = fun(N) giving the deterministic number of iterations for the first pass. |

`alpha` |
when using default |

`compute_confidence` |
if |

`verbose` |
to control level of print (either 0, or 1, or 2). |

`chi2` |
for a chi2 test on the number of events. |

`breaks` |
for the final histogram is |

`...` |
further arguments past to |

`MP`

is a wrap up of `IRW`

for probability and quantile estimation.
By construction, the several calls to `IRW`

are parallel (foreach)
and so is the algorithm. Especially, with `N.batch=1`

, this is the Last Particle
Algorithm, which is a specific version of `SubsetSimulation`

with `p_0 = 1-1/N`

.
However, note that this algorithm not only gives a quantile or a probability estimate
but also an estimate of the whole cdf until the given threshold `q`

.

The probability estimator only requires to generate several random walks as it is the estimation of the parameter of a Poisson random variable. The quantile estimator is a little bit more complicated and requires a 2-passes algorithm. It is thus not exactly fully parallel as cluster/cores have to communicate after the first pass. During the first pass, particles are moved a given number of times, during the second pass particles are moved until the farthest event reach during the first pass. Hence, the random process is completely simulated until this given state.

For an easy user experiment, all the parameters are defined by default with the optimised values
as described in the reference paper (see References below) and a typical use will only specify
`N`

and `N.batch`

.

An object of class `list`

containing the outputs described below:

`p` |
the estimated probability or the reference for the quantile estimate. |

`q` |
the estimated quantile or the reference for the probability estimate. |

`ecdf_MP` |
the empirical cdf. |

`L_max` |
the farthest state reached by the random process. Validity range
for the |

`times` |
the |

`Ncall` |
the total number of calls to the |

`particles` |
the |

`LSF_particles` |
the value of the |

`moves` |
a vector containing the number of moves for each batch. |

`p_int` |
a 95% confidence intervalle on the probability estimate. |

`q_int` |
a 95% confidence intervall on the quantile estimate. |

`chi2` |
the output of the chisq.test function. |

The `alpha`

parameter is set to 0.05 by default. Indeed it should not be
set too small as it is defined approximating the Poisson distribution with the Gaussian one.
However if no estimate is produce then the algorithm can be restarted for the few missing events.
In any cases, setting `Niter_1fold = -N/N.batch*log(p)`

gives 100% chances to produces
a quantile estimator.

Clement WALTER clement.walter@cea.fr

A. Guyader, N. Hengartner and E. Matzner-Lober:

*Simulation and estimation of extreme quantiles and extreme probabilities*

Applied Mathematics \& 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.

E. Simonnet:

*Combinatorial analysis of the adaptive last particle method*

Statistics and Computing, 1-20.

`SubsetSimulation`

`MonteCarlo`

`IRW`

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 | ```
## Not run:
# Estimate some probability and quantile with the parabolic lsf
p.est <- MP(2, kiureghian, N = 100, q = 0) # estimate P(lsf(X) < 0)
p.est <- MP(2, kiureghian, N = 100, q = 7.8, lower.tail = FALSE) # estimate P(lsf(X) > 7.8)
q.est <- MP(2, kiureghian, N = 100, p = 1e-3) # estimate q such that P(lsf(X) < q) = 1e-3
q.est <- MP(2, kiureghian, N = 100, p = 1e-3, lower.tail = FALSE) # estimate q such
# that P(lsf(X) > q) = 1e-3
# plot the empirical cdf
plot(xplot <- seq(-3, p.est$L_max, l = 100), sapply(xplot, p.est$ecdf_MP))
# check validity range
p.est$ecdf_MP(p.est$L_max - 1)
# this example will fail because the quantile is greater than the limit
tryCatch({
p.est$ecdf_MP(p.est$L_max + 0.1)},
error = function(cond) message(cond))
# Run in parallel
library(doParallel)
registerDoParallel()
p.est <- MP(2, kiureghian, N = 100, q = 0, N.batch = getDoParWorkers())
## End(Not run)
``` |

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.

All documentation is copyright its authors; we didn't write any of that.