IRW: Increasing Randow Walk

In clemlaflemme/mistral: Methods in Structural Reliability

Description

Simulate the increasing random walk associated with a real-valued continuous random variable.

Usage

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

Arguments

dimension	dimension of the input space.
lsf	limit state function.
N	number of particules.
q	level until which the randow walk is to be generated.
Nevent	the number of desired events.
X	to start with some given particles.
y	value of the lsf on X.
K	kernel transition for conditional generations.
burnin	burnin parameter.
sigma	radius parameter for K.
last.return	if the last event should be returned.
use.potential	tu use a ‘potential’ matrix to select starting point not directly related to the sample to be moved with the MH algorithm.
plot	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.
plot.lsf	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.
print_plot	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.
output_dir	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.
plot.lab	the x and y labels for the plot

Details

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 μ^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.

Value

An object of class list containing the following data:

L	the events of the random walk.
M	the total number of iterations.
Ncall	the total number of calls to the lsf.
X	a matrix containing the final particles.
y	the value of lsf on X.
q	the threshold considered when generating the random walk.
Nevent	the target number of events when generating the random walk.
Nwmoves	the number of rejected transitions, ie when the proposed point was not stricly greater/lower than the current state.
acceptance	a vector containing the acceptance rate for each use of the MH algorithm.

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; see examples in MonteCarlo.

Author(s)

Clement WALTER [e-mail]

References

• 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 \& 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 \& Optimization, 64(2), 171-196.

See Also

MP

Examples

# 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)
## Not run: 
fail.samp <- IRW(2, lsf, q = 0, N = 10, plot = TRUE)

## End(Not run)

clemlaflemme/mistral documentation built on Jan. 3, 2020, 9:13 a.m.