MC_IS: Monte Carlo Simulation with Importance Sampling

In TesiproV: Reliability Analysis Methods for Structural Engineering

Monte Carlo Simulation with Importance Sampling

Description

Performs Monte Carlo simulation with importance sampling (IS) to estimate the failure probability of a single limit-state function or of a system of limit-state functions.

Failure is defined as g(x) < 0. Samples with g(x) = 0 are counted as safe.

The samples are drawn from a sampling density h in standard normal space that is centred at FORM design points. The sampling is corrected by the likelihood ratio w(u) = \varphi(u) / h(u), where \varphi is the standard normal density.

Usage

MC_IS(
  lsf,
  lDistr,
  cov_user = 0.05,
  n_batch = 5000,
  n_max = 1e+06,
  use_threads = NULL,
  backend = NULL,
  sys_type = "serial",
  dataRecord = TRUE,
  beta_l = 100,
  densityType = "norm",
  dps = NULL,
  debug.level = 0,
  seed = NULL,
  adaptive_alpha = FALSE,
  alpha_update_rate = 0.1,
  adaptive_batch = FALSE,
  batch_control = list(),
  min_adapt_samples = 10000,
  alpha_min = 0.02,
  stability_mode = c("robust", "fast")
)

Arguments

lsf	A limit-state function or a list of limit-state functions for a system. Must follow the LSF interface described below.
lDistr	List of marginal distribution objects (funlists with $d, $p, $q, $r). For a system, a list with one such list per limit-state function.
cov_user	Target coefficient of variation of the estimator.
n_batch	Number of samples per batch.
n_max	Maximum total number of samples (see stopping criterion).
use_threads	Number of chunks per batch. With backend = "parallel" also the number of cores or workers.
backend	Parallel backend, "future" or "parallel". If NULL (default), "parallel" with one worker is used, i.e. the calculation runs sequentially on every system.
sys_type	"serial" or "parallel". Only used for systems.
dataRecord	Logical. If TRUE, the results of each batch are returned in data. Single limit-state function only.
beta_l	Currently without effect.
densityType	Currently without effect. A normal sampling density in standard normal space is always used.
dps	Optional design points in physical space. If supplied, FORM is skipped for the components. Single limit-state function: a numeric vector with one value per variable (or one value for all). System: a list with one such vector per limit-state function. A numeric vector with one value per limit-state function is accepted if each function has exactly one variable.
debug.level	Verbosity level: 0 = silent, 1 = summary, 2 = detailed output.
seed	Optional integer seed for reproducible results.
adaptive_alpha	Logical. Adaptive mixture weights are not implemented. Without effect for a single limit-state function. For serial multimodal systems the calculation stops with an error once min_adapt_samples is exceeded.
alpha_update_rate	Reserved for adaptive_alpha, currently without effect.
adaptive_batch	Logical. Systems only. If TRUE, the batch size is adapted from the current effective number of failure samples, so that cov_user is expected to be reached after K_future further batches.
batch_control	List for adaptive_batch with the elements n_min and n_max (bounds of the batch size) and K_future (number of remaining batches). Systems only.
min_adapt_samples	Reserved for adaptive_alpha (see there).
alpha_min	Reserved for adaptive_alpha, currently without effect.
stability_mode	"robust" or "fast". Without influence on the result, all sums are accumulated in the log domain. For serial multimodal systems "fast" additionally returns the effective sample size per mixture component (ESS_modes).

Details

Estimator

The failure probability and the variance of the estimator are

\hat{P}_f = \frac{1}{N} \sum_{i=1}^{N} I_i w_i, \qquad \widehat{\mathrm{Var}}(\hat{P}_f) = \frac{1}{N}\left(\frac{1}{N}\sum_{i=1}^{N} I_i w_i^2 - \hat{P}_f^2\right),

where I_i = 1 if sample i lies in the failure domain and I_i = 0 otherwise. The coefficient of variation is \sqrt{\widehat{\mathrm{Var}}(\hat{P}_f)} / \hat{P}_f.

Single limit-state function

The failure domain is \{g(x) < 0\}. The sampling density is a multivariate normal density centred at the FORM design point u^* with the covariance matrix

\Sigma = I + (\beta^2 - 1)\, \alpha \alpha^T, \qquad \alpha = u^* / \|u^*\|,

where \beta is the FORM reliability index, or \|u^*\| if it is not available. The variance is \beta^2 in the direction of \alpha and 1 perpendicular to it. FORM is only used to locate the centre. A design point from a FORM run that did not converge is accepted as long as it is finite.

Systems

For limit-state functions g_1, \dots, g_m the failure domain is \bigcup_k \{g_k(x) < 0\} for a serial system and \bigcap_k \{g_k(x) < 0\} for a parallel system. Basic variables that appear in several limit-state functions are identified by their name and enter the simulation only once.

The sampling density is a mixture of unit-variance normal densities

h(u) = \sum_{j} \lambda_j \, \varphi(u - u_j^*), \qquad \sum_j \lambda_j = 1.

• Serial systems: the centres u_j^* are the FORM design points of the components. The weights are proportional to \Phi(-\beta_j) if all FORM results are finite, otherwise they are equal. If the largest distance between the design points is smaller than 1.5 and the smallest cosine between their gradient directions is larger than 0.85, a single normal density at the mean of the design points is used instead of the mixture.

• Parallel systems: FORM is run for the joint failure domain \max_k g_k(x) < 0. If its design point is available, it is a centre with weight 0.7 and the remaining weight is divided equally among the component design points. Otherwise the component design points are used with equal weights.

Stopping criterion

For a single limit-state function the simulation stops when the coefficient of variation falls below cov_user or when n_max samples are reached. The last batch is reduced so that the total does not exceed n_max.

For systems the simulation stops when the coefficient of variation falls below cov_user, the effective number of failure samples exceeds 30 (serial) or 50 (parallel) and more than 5 * n_batch samples have been used. It also stops when n_max is reached. The batches are not reduced, so n_mc can exceed n_max by less than one batch. For parallel systems a warning is issued if the effective number of failure samples is below 50 at the end.

Numerical stability

The weighted sums are accumulated in the log domain. The result includes log_pf. For extremely small probabilities pf may underflow to zero, beta is then computed from log_pf and can remain finite. This applies to single limit-state functions and to systems.

Automatic detection of vectorized limit-state functions

Before the simulation the limit-state function is called once with a matrix of two rows (the design point in physical space). If it returns a finite numeric vector of length 2, it is evaluated for whole batches at once. Otherwise it is evaluated row by row with apply(). Scalar limit-state functions therefore need not be changed.

Parallel execution

The samples of each batch are generated in the master process and split into use_threads chunks, which are evaluated in parallel. The workers do not generate random numbers.

• backend = "future": the chunks are evaluated with future.apply::future_lapply() under the active future plan. The function never changes the plan. If the plan is sequential and use_threads > 1, a message is shown once per session and the calculation runs sequentially. The number of workers is determined by the plan. use_threads only determines the number of chunks.

• backend = "parallel": on Unix-like systems parallel::mclapply() with use_threads cores is used. On Windows a PSOCK cluster with use_threads workers is created for the call and stopped afterwards. With use_threads = 1 lapply() is used.

Reproducibility

With a fixed seed the random numbers are generated from a L'Ecuyer-CMRG stream in the master process. The sample therefore does not depend on use_threads. Results for different values of use_threads agree up to rounding errors, because the sums are accumulated per chunk. The function sets RNGkind("L'Ecuyer-CMRG") and changes .Random.seed in the global environment.

Value

For a single limit-state function a list with

• method - "MCIS_Single"

• beta - reliability index -\Phi^{-1}(\hat{P}_f)

• pf - estimated failure probability

• log_pf - natural logarithm of pf

• log_var - natural logarithm of var

• FORM_beta, FORM_pf - FORM results (NA if dps was supplied)

• design_points - centre of the sampling density in standard normal space

• var - variance of the estimator

• cov_mc - coefficient of variation of the estimator

• cov_user - target coefficient of variation

• ESS - effective sample size (\sum w_i)^2 / \sum w_i^2

• ESS_ratio - ESS / n_mc

• n_mc, n_max - used and maximum number of samples

• data - results per batch (if dataRecord = TRUE)

• runtime - runtime as returned by proc.time()

For a system a list with

• method - "MCIS_System"

• sys_type - "serial" or "parallel"

• beta, pf, var, cov_mc - as above

• log_pf, log_var - as above

• RSE_est - identical to cov_mc

• ESS - effective number of failure samples (\sum I_i w_i)^2 / \sum (I_i w_i)^2

• ESS_modes - effective sample size per mixture component (serial systems with stability_mode = "fast", otherwise NULL)

• n_mc - number of samples used

Limit-State Function (LSF) Interface

The limit-state function follows the sign convention g(x) > 0 safe, g(x) = 0 limit state and g(x) < 0 failure. It must use one of the following signatures:

function(x)

A single numeric vector x containing all basic variables in the order defined in lDistr. This form also allows vectorized evaluation and is recommended.

function(Z, Fy, M, ...)

Explicitly named scalar arguments corresponding to the variable names defined in the probabilistic model.

The use of function(...) as the sole argument is not supported, because it leads to inconsistent behaviour between FORM-based and Monte Carlo-based reliability methods.

Windows users

On Windows, backend = "future" with future::plan(multisession, workers = ...) is recommended. backend = "parallel" creates a new PSOCK cluster for every call and issues a warning.

Author(s)

(C) 2021-2026 K. Nille-Hauf, T. Feiri, M. Ricker, T. Lux – Hochschule Biberach (until 2022), TU Dortmund University - Chair of Structural Concrete (since 2023)

References

Ditlevsen, O., & Madsen, H. O. (1996). Structural Reliability Methods. Wiley.

Melchers, R. E., & Beck, A. T. (2018). Structural Reliability Analysis and Prediction. Wiley.

Spaethe, G. (1991). Die Sicherheit tragender Baukonstruktionen. Springer.

TesiproV documentation built on Oct. 28, 2026, 5:13 p.m.