2. Monte-Carlo Sampling and Importance Sampling

#math/distribution-transform#math/Importance-sampling#math/monte-carlo-sampling

1. Basic Concepts

(1) Sampling Preliminaries

1) Distribution transformation

In the sampling functions, we call $f_X(x)$ as the probability density. In most of cases, we firstly sample from the range (0, 1) and then transform it by some distribution function.

To sample in numpy, we can use rng = np.random.default_rng and then use rng.uniform to sample a uniform distribution. We set

$$\begin{array}{}\text{(1.1.1)}& U\sim \mathcal{U}(0,1)U(x)=y\end{array}$$

Then we use an exponential distribution with parameter $\lambda$, the distribution function is :

$$\begin{array}{}\text{(1.1.2)}& f_X=\lambda e^{-\lambda x}\end{array}$$

the cumulative distribution function would be :

$$\begin{array}{}\text{(1.1.3)}& F_X(x)=z=1-e^{-\lambda x}{F}^{-1}(x)=-\frac{\ln (1-F_X(x))}{\lambda }\end{array}$$

in order to transform a $U$ into $F$, we need to find a transform $K$ to map $x_1$ is mapped to $x_2$, i.e.,

$$\begin{array}{}\text{(1.1.4)}& U(x_1)=F(x_2)\end{array}$$

Since $x_1=U^{-1}(y)$ Then the transforming is :

$$\begin{array}{}\text{(1.1.5)}& \boxed{{\displaystyle x_2=F^{-1}(U(x_1))}}\end{array}$$

here $U$ is cumulative distribution function of original variable, $F^{-1}$ is cumulative distribution of target variable.

As a general case, the following part gives the transform from a normal distribution to a exponential distribution :

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

# Number of samples
n_samples = 10_000

# Create a random number generator
rng = np.random.default_rng(seed=42)

# Step 1: Sample U from the uniform distribution (0, 1)
u = rng.normal(0.0, 1.0, size=n_samples)

# Step 2: Transform U into an exponential random variable
lam = 1.0  # Exponential distribution rate parameter
x = -np.log(1.0 - norm.cdf(u)) / lam

# Step 3: Plot the samples as a bar-style histogram
x_max = np.percentile(x, 99)  # Avoid displaying extreme outliers
bins = np.linspace(0, x_max, 40)
# Theoretical exponential probability density function
x_theory = np.linspace(0, x_max, 500)
pdf = lam * np.exp(-lam * x_theory) 

plt.hist(
    x,
    bins=bins,
    density=True,
    alpha=0.75,
    color="steelblue",
    edgecolor="black",
    label="Generated samples"
)
plt.plot(
    x_theory,
    pdf,
    color="red",
    linewidth=2,
    label=r"Theoretical PDF: $\lambda e^{-\lambda x}$"
)

plt.xlabel("x")
plt.ylabel("Probability density")
plt.title("Sampling an Exponential Distribution by Inverse CDF")
plt.legend()
plt.tight_layout()
plt.show()

Which can yield the distribution graph as follows :

2) Reliability Analysis & Monte Carlo sampling

Firstly, to define whether a structure is safe or failure, we use a performance function :

$$\begin{array}{}\text{(1.1.6)}& g(x)=\{\begin{array}{l}\ge 0\text{safe}\\ <0\text{failure}\end{array}\end{array}$$

And we define the failure region as $F$, The Monte Carlo sampling is sampling directly by $f_X(x)$, then the failure probability is :

$$\begin{array}{}\text{(1.1.7)}& P_f={\int }_F{f}_X(x)dx={\int }_{R^n}{I}_{F}f(x)dx{I}_f=\{\begin{array}{l}1\text{failure}\\ 0\text{safe}\end{array}\end{array}$$

where $I_F$ is the indicator function of failure.

We note since we already use $f(x)$ as the sampling density, the failure probability is defined by :

$$\begin{array}{}\text{(1.1.8)}& {\hat{P}}_f=E({\hat{P}}_f)=\frac{1}{N}\sum _{j=1}^n{I}_F(x_j)=\frac{N_F}{N}{x}_j\sim f(x)\end{array}$$

The variance of the failure probability can be calculated as :

$$\begin{array}{}\text{(1.1.9)}& V[{\hat{P}}_f]=V[\frac{1}{N}\sum _{i=1}^n{I}_F(x_i)]\approx \frac{({\hat{P}}_f-{\hat{P}}_f^2)}{N-1}\end{array}$$

The coefficient of variation (CoV) of the Monte Carlo estimator is :

$$\begin{array}{}\text{(1.1.10)}& \mathrm{CoV}[{\hat{P}}_f]=\frac{\sqrt{V[{\hat{P}}_f]}}{{\hat{P}}_f}\approx \sqrt{\frac{1-{\hat{P}}_f}{(N-1){\hat{P}}_f}}.\end{array}$$

(2) Importance sampling

1) Computation of Failure probability

We know in the Monte-Carlo sampling method, we define the performance function $g(X)$ (denoting the structural fidelity),

$$\begin{array}{}\text{(1.2.1)}& P_f={\int }_F{I}_F{f}_X(x)dx\underset{\text{for independent variables}}{\underbrace{f_X(x)=\prod _{i=1}^n{f}_{X_i}(x_i)}}\end{array}$$

where $f_X(x)$ is the unite sampling density distribution function. Which is the variable distribution in the sampling case (In engineering, all variables often satisfies the normal distribution). We also note $i$ is for the variable, $j$ is for the sampling identity :

$$\begin{array}{}\text{(1.2.2)}& f_X(x_j)=\{\begin{array}{l}1\text{ failure }\\ 0\text{ valid}\end{array}\end{array}$$

However, we can still use another sampling function $h_X(x)$, then the failure probability is given by :

$$\begin{array}{}\text{(1.2.3)}& {\hat{P}}_f={\int }_{R^n}{I}_F(x)\frac{f_X(x)}{h_X(x)}{h}_X(x)dx=\frac{1}{N}\sum _{j=1}^N[I_F(x_j)\frac{f_X(x_j)}{h_X(x_j)}],x_j\sim h_X.\end{array}$$

This is called Importace sampling^[1][2]

2) Mean, standard deviation, and coefficient of variation

Define the importance weight and the contribution of sample $j$ as (weighted failure probability):

$$\begin{array}{}\text{(1.2.4)}& w(x)=\frac{f_X(x)}{h_X(x)},Y_j=I_F(x_j)w(x_j),x_j\sim h_X.\end{array}$$

The importance-sampling estimator is ${\hat{P}}_f=N^{-1}\sum _{j=1}^N{Y}_j$. Provided $h_X(x)>0$ wherever $I_F(x)f_X(x)>0$, its mean is

$$\begin{array}{}\text{(1.2.5)}& E_h[{\hat{P}}_f]=E_h[Y]={\int }_{R^n}{I}_F(x)\frac{f_X(x)}{h_X(x)}{h}_X(x)dx=P_f.\end{array}$$

Thus, the estimator is unbiased. Its variance (theoretical variance) is :

$$\begin{array}{}\text{(1.2.6)}& {\mathrm{Var}}_h[{\hat{P}}_f]=\frac{{\mathrm{Var}}_h[Y]}{N}=\frac{1}{N}[{\int }_{R^n}{I}_F(x)\frac{f_X^2(x)}{h_X(x)}dx-P_f^2].\end{array}$$

Therefore, the standard deviation (standard error) is :

$$\begin{array}{}\text{(1.2.7)}& {\mathrm{Std}}_h[{\hat{P}}_f]=\sqrt{\frac{E_h[Y^2]-P_f^2}{N}}.\end{array}$$

For a simulation estimate, let $\overline{Y}={\hat{P}}_f$. An unbiased estimator of the variance and the corresponding estimated standard deviation (sample-based estimator) are defined as :

$$\begin{array}{}\text{(1.2.8)}& \widehat{\mathrm{Var}}[{\hat{P}}_f]=\frac{1}{N(N-1)}\sum _{j=1}^N(Y_j-\overline{Y}{)}^2,\widehat{\mathrm{Std}}[{\hat{P}}_f]=\sqrt{\widehat{\mathrm{Var}}[{\hat{P}}_f]}.\end{array}$$

We also note in (1.2.8), the term of $\widehat{Var}$ can be computed by $x$ as following equation :

$$\begin{array}{}\text{(1.2.9)}& \widehat{\mathrm{Var}}[{\hat{P}}_f]=\frac{1}{N(N-1)}\sum _{j=1}^N{(I_f\frac{f_X}{h_X}-{\hat{P}}_f)}^2=\frac{1}{N(N-1)}[\sum _{j=1}^N{I}_f\frac{f_X^2}{h_X^2}-{\hat{P}}_f\underset{N{\hat{P}}_f}{\underbrace{\sum _{j=1}^N{I}_f\frac{f_X}{h_X}}}+N{\hat{P}}_f^2]=\frac{1}{N-1}[\frac{1}{N}\sum _{j=1}^N{I}_f\frac{f_X^2}{h_X^2}-{\hat{P}}_f^2]\end{array}$$

Finally, for $P_f>0$, the coefficient of variation is :

$$\begin{array}{}\text{(1.2.10)}& \mathrm{CoV}[{\hat{P}}_f]=\frac{{\mathrm{Std}}_h[{\hat{P}}_f]}{P_f}=\sqrt{\frac{E_h[Y^2]-P_f^2}{N{P}_f^2}},\end{array}$$

and its practical estimate is

$$\begin{array}{}\text{(1.2.11)}& \widehat{\mathrm{CoV}}[{\hat{P}}_f]=\frac{\widehat{\mathrm{Std}}[{\hat{P}}_f]}{{\hat{P}}_f}.\end{array}$$

A good importance density $h_X$ assigns substantial probability to the failure region and makes $Y=I_F{f}_X/h_X$ nearly constant, thereby reducing the standard deviation and CoV.

3) Optimal Importance sampling function

In the structural reliability analysis, we often sample near the failure point, we define as the "design point", which is $x_0$ here. We want to sample more failure points near the failure surface, so the more it easy to fail, we are more willing to sample it. The optimal importance sampling function can be expressed by $I_F$, as :

$$\begin{array}{}\text{(1.2.12)}& h_X^{opt}(x)=\frac{I_F(x)f_X(x)}{{\hat{P}}_f}\end{array}$$

In practical case, we often use the normal distribution near the design point, i.e.,

$$\begin{array}{}\text{(1.2.13)}& \boxed{{\displaystyle h_X(x)=\frac{1}{\sqrt{2\pi }{\sigma }_s}\exp {\left(\frac{x-x_0}{{\sigma }_s}\right)}^2}}\end{array}$$

Reference

---

1. https://en.wikipedia.org/wiki/Importance_sampling ↩︎

2. 吕震宙, 宋述芳, 李璐伟等, 结构/机构可靠性设计基础 [M] 西北工业大学出版社, 2019.3 ↩︎