Simulate the increasing random walk associated with a real-valued continuous random variable.
IRW(
dimension,
lsf,
N = 10,
q = Inf,
Nevent = Inf,
X,
y = lsf(X),
K,
burnin = 20,
sigma = 0.3,
last.return = TRUE,
use.potential = TRUE,
plot = FALSE,
plot.lsf = FALSE,
print_plot = FALSE,
output_dir = NULL,
plot.lab = c("x_1", "x_2")
)
An object of class list
containing the following data:
the events of the random walk.
the total number of iterations.
the total number of calls to the lsf
.
a matrix containing the final particles.
the value of lsf
on X
.
the threshold considered when generating the random walk.
the target number of events when generating the random walk.
the number of rejected transitions, ie when the proposed point was not stricly greater/lower than the current state.
a vector containing the acceptance rate for each use of the MH algorithm.
dimension of the input space.
limit state function.
number of particules.
level until which the randow walk is to be generated.
the number of desired events.
to start with some given particles.
value of the lsf
on X.
kernel transition for conditional generations.
burnin parameter.
radius parameter for K
.
if the last event should be returned.
tu use a ‘potential’ matrix to select starting point not directly related to the sample to be moved with the MH algorithm.
if TRUE
, the algorithm plots the evolution of the particles. This
requieres to evaluate the lsf
on a grid and is only for visual purpose.
a boolean indicating if the lsf
should be added to the
plot. This requires the evaluation of the lsf
over a grid and
consequently should be used only for illustation purposes.
if TRUE, print the updated plot after each iteration. This might
be slow; use with a small N
. Otherwise it only prints the final plot.
if plots are to be saved in pdf in a given directory. This will
be pasted with ‘_IRW.pdf’. Together with print_plot==TRUE
this will
produce a pdf with a plot at each iteration, enabling ‘video’ reconstitution
of the algorithm.
the x and y labels for the plot
Clement WALTER clementwalter@icloud.com
This function lets generate the increasing random walk associated with a continous
real-valued random variable of the form Y = lsf(X)
where X
is
vectorial random variable.
This random walk can be associated with a Poisson process with parameter
N
and hence the number of iterations before a given threshold q
is directly related to P[ lsf(X) > q]. It is the core tool of algorithms
such as nested sampling, Last Particle Algorithm or Tootsie Pop Algorithm.
Bascially for N = 1
, it generates a sample \(Y = lsf(X)\) and iteratively
regenerates greater than the found value: \(Y_{n+1} \sim \mu^Y( \cdot \mid Y > Y_n\). This
regeneration step is done with a Metropolis-Hastings algorithm and that is why it is usefull
to consider generating several chains all together (N > 1
).
The algorithm stops when it has simulated the required number of events Nevent
or when
it has reached the sought threshold q
.
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.
C. Walter:
Point Process-based Monte Carlo estimation
Statistics and Computing, in press, 1-18.
arXiv preprint arXiv:1412.6368.
J. Skilling:
Nested sampling for general Bayesian computation
Bayesian Analysis, 1(4), 833-859.
M. Huber and S. Schott:
Using TPA for Bayesian inference
Bayesian Statistics 9, 9, 257.
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.
MP
# Get faililng samples for the kiureghian limit state function
# Failure is defined as lsf(X) < 0 so we have to invert the lsf
lsf <- function(x) -1*kiureghian(x)
if (FALSE) {
fail.samp <- IRW(2, lsf, q = 0, N = 10, plot = TRUE)
}
Run the code above in your browser using DataLab