Probability transformations and copulas

Probability transformation

Classical FORM uses independent standard normal coordinates. A probability transformation maps the joint law of the physical variables into that space. Generalized Nataf also permits spherical non-normal standard spaces, provided the reliability calculation uses the corresponding probability law.

Transformation of dependent random variables using Nataf approach

One method to handle this is using the Nataf joint distribution model, if the marginal cdfs are known. [Baker2010] The correlated random variables \(\mathbf{X} = ( X_1 , \dots , X_n )\) with the correlation matrix \(\mathbf{R}\) can be transformed by

(1)\[ y_i=\Phi^{-1}\left(F_{X_{i}}(x_i)\right) \qquad i = 1,\dots,n\]

into normally distributed random variables \(\mathbf{Y}\) with zero means and unit variance, but still correlated with \(\mathbf{R}_0\) . Nataf’s distribution for \(\mathbf{X}\) is obtained by assuming that \(\mathbf{Y}\) is jointly normal. [Liu1986]

The correlation coefficients for \(\mathbf{X}\) and \(\mathbf{Y}\) are related by

(2)\[ \rho_{X_i,X_j} = \int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty} \left(\frac{x_i-\mu_{X_i}}{\sigma_{X_i}}\right) \left(\frac{x_j-\mu_{X_j}}{\sigma_{X_j}}\right) \frac{1}{2\pi \sqrt{1-\rho_{Y_i,Y_j}^2}} \exp\left(-\frac{y_i^2-2\rho_{Y_i,Y_j}y_iy_j+y_j^2}{2(1-\rho_{Y_i,Y_j}^2)}\right) d y_i d y_j\]

Once this is done, the transformation from the correlated normal random variables \(\mathbf{Y}\) to uncorrelated normal random variables \(\mathbf{Z}\) is addressed. Hence, the transformation is

(3)\[ \mathbf{z}=\mathbf{L}_0^{-1}\mathbf{y} \quad \Leftrightarrow \quad \mathbf{y} = \mathbf{L}_0\mathbf{z}\]

where \(\mathbf{L}_0\mathbf{L}_0^T=\mathbf{R}_0\) is the Cholesky factorization of the correlation matrix of \(\mathbf{Y}\). The Jacobian matrix, denoted by \(\mathbf{J}\), for the transformation is given by

(4)\[ \mathbf{J}_{ZX} = \frac{\partial \mathbf{z}}{\partial \mathbf{x}} = \mathbf{L}_0^{-1}\text{diag} \left(\frac{f_{X_i}(x_i)}{\varphi(y_i)}\right)\]

This approach is useful when the marginal distribution for the random variables \(\mathbf{X}\) is known and the knowledge about the variables dependence is limited to correlation coefficients. [Baker2010] [DerKiureghian2006]

Transformation of dependent random variables using Rosenblatt approach

An alternative to the Nataf approach is to consider the joint pdf of \(\mathbf{X}\) as a product of conditional pdfs.

(5)\[ f_{\mathbf{X}}(\mathbf{x}) = f_{X_1}(x_1) f_{X_2|X_1}(x_2|x_1) \dots f_{X_n|X_1,\dots,X_{n-1}}(x_n|x_1,\dots,x_{n-1})\]

As a result of the sequential conditioning in the pdf, the conditional cdfs are given for \(i \in [1,n]\)

(6)\[ F_{X_i|X_1,\dots,X_{i-1}}(x_i|x_1,\dots,x_{i-1}) = \int_{-\infty}^{x_i} f_{X_i|X_1,\dots,X_{i-1}}(x_i|x_1,\dots,x_{i-1}) d x_i\]

These conditional distributions for the random variables \(\mathbf{X}\) can be transformed into standard normal marginal distributions for the variables \(\mathbf{Z}\), using the so called Rosenblatt transformation [Rosenblatt1952], suggested by Hohenbichler and Rackwitz (1981).

(7)\[\begin{split} \begin{split} z_1 &= \Phi^{-1}\left( F_{X_1}(x_1) \right)\\ z_2 &= \Phi^{-1}\left( F_{X_2|X_1}(x_2|x_1) \right)\\ &\vdots\\ z_n &= \Phi^{-1}\left( F_{X_n|X_1,\dots,X_{n-1}}(x_n|x_1,\dots,x_{n-1}) \right) \end{split}\end{split}\]

The Jacobian of this transformation is a lower triangular matrix having the elements [Baker2010]

(8)\[\begin{split} \left[\mathbf{J}_{ZX}\right]_{i,j} = \frac{\partial z_i}{\partial x_j} = \begin{cases}\displaystyle \frac{1}{\varphi(z_i)}\frac{\partial}{\partial x_j} F_{X_i|X_1,\dots,X_{i-1}}(x_i|x_1,\dots,x_{i-1}) & i \geq j\\ 0 & i < j \end{cases}\end{split}\]

Here \(\varphi\) is the standard-normal density, and the triangular structure follows the chosen conditioning order.

In some cases the Rosenblatt transformation cannot be applied, because the required conditional pdfs cannot be provided. In this case other transformations may be useful, for example Nataf transformation. [Faber2009]

Copulas as joint distribution specifications

For continuous marginals, dependence is specified by a copula \(C\):

\[F_{\mathbf X}(\mathbf x)=C(F_1(x_1),\ldots,F_n(x_n)),\qquad f_{\mathbf X}(\mathbf x)=c(F_1(x_1),\ldots,F_n(x_n))\prod_i f_i(x_i).\]

The copula and marginals define the joint law before a transformation is chosen. Physical Pearson correlations alone do not determine that law. PySTRA’s legacy Pearson interface assumes a Gaussian copula and calibrates its latent correlation using (2). An explicit GaussianCopula already specifies the latent correlation, so that calibration is bypassed. Changing to Student-t or Frank changes the dependence law, rather than just changing a numerical factorization.

For Gaussian and Student-t copulas, the pairwise Kendall rank correlations and latent unit-diagonal shape entries satisfy

\[\tau_{ij}=\frac{2}{\pi}\arcsin R_{ij}.\]

The resulting full shape matrix must be positive definite. Arbitrary physical marginals can be used with the same copula. See Dependence and probability transformations for sampling and density examples, and Copulas and isoprobabilistic transformations for the supported specification interface.

Generalized Nataf for elliptical copulas

Following [LebrunDutfoy2009a], let \(H\) be the univariate CDF of a centered elliptical representative with unit-diagonal shape matrix \(R=LL^T\). Generalized Nataf is

\[w_i=H^{-1}(F_i(x_i)),\qquad \mathbf v=L^{-1}\mathbf w, \qquad x_i=F_i^{-1}(H((L\mathbf v)_i)).\]

The transformed distribution is spherical. Its coordinates are independent only in the Gaussian case within these Gaussian/Student-t families. For Student-t, \(H=T_\nu\) and \(\mathbf V\) is spherical Student-t with identity shape. When \(\nu>2\), its covariance is \(\nu/(\nu-2)I\); for smaller degrees of freedom, covariance does not exist. Identity shape therefore must not be interpreted as independence. The Nataf Jacobian is

\[J_{VX}=L^{-1}\operatorname{diag}\left(\frac{f_i(x_i)}{h(w_i)}\right),\]

where \(h=H'\). Spherical symmetry makes every unit projection have CDF \(H\). Hence a tangent failure half-space at signed geometric distance \(\beta\) has probability

\[P_{f,\mathrm{FORM}}=H(-\beta).\]

For Student-t, PySTRA uses \(T_\nu(-\beta)\), while FORMResult.beta reports \(-\Phi^{-1}(P_f)\). SORM, system FORM, simulation methods and the current Strong Maximum Test require independent normal coordinates and reject spherical Student-t Nataf space. A normal Rosenblatt mapping supports these analyses with the same Student-t copula.

Rosenblatt conditioning and approximation invariance

A copula supplies the conditional CDFs required by Rosenblatt. For an order \(\pi\), define \(p_i=F_i(x_i)\) and

\[q_{\pi_k}=C_{\pi_k\mid\pi_1,\ldots,\pi_{k-1}} (p_{\pi_k}\mid p_{\pi_1},\ldots,p_{\pi_{k-1}}),\qquad u_{\pi_k}=\Phi^{-1}(q_{\pi_k}).\]

The \(q_i\) are independent uniforms and the \(u_i\) are independent standard normals. PySTRA keeps coordinate arrays in original marginal order; u[order[k]] denotes the k-th conditional innovation.

For a Student-t copula, reorder the latent vector and shape first, and write \(\mathbf z=L^{-1}\mathbf w\). With one-based conditional index \(k\),

\[s_k^2=\frac{\nu+\sum_{j<k}z_j^2}{\nu+k-1},\qquad q_{\pi_k}=T_{\nu+k-1}(z_k/s_k).\]

The inverse constructs each \(z_k=s_kT_{\nu+k-1}^{-1}(q_{\pi_k})\) sequentially, then recovers the latent vector and physical marginals. This additional conditional scaling removes the dependence that remains after Student-t Nataf whitening.

For a Gaussian copula, Cholesky Nataf and Rosenblatt coincide for a fixed order. Different orders are related by an orthogonal change of normal coordinates; the optimum FORM distance and probability are invariant. Numerical optimization can still converge to different local points. For non-Gaussian copulas, transformations need not be related orthogonally: exact probability remains invariant, but FORM’s tangent approximation can change with conditioning order [LebrunDutfoy2009b]. All system components must use the same mapping and order before their normal directions are compared.

The bivariate Frank copula used in that paper is

\[C_\theta(p_1,p_2)=-\frac{1}{\theta}\log\left[ 1+\frac{(e^{-\theta p_1}-1)(e^{-\theta p_2}-1)}{e^{-\theta}-1}\right],\]

with the independence limit at \(\theta=0\). Its conditional CDF is \(C_{2\mid1}=\partial C_\theta/\partial p_1\). For exponential rates 1 and 3, \(\theta=10\), and failure event \(8X_1+2X_2-1\leq0\), the tutorial reproduces FORM probabilities about 0.107 and 0.122 for the two orders, compared with direct integration of about 0.1038. This difference is approximation error, not a different underlying probability for each order.

Very small probabilities

Every marginal mapping above passes through a probability, \(\Phi(u)\) or \(F_X(x)\). In double precision \(\Phi(u)\) rounds to one above \(u \approx 8.3\), so \(F_X^{-1}(\Phi(u))\) returns the upper end of the support, which is infinite for most load models. Line sampling and importance sampling at high reliability visit such points routinely, and FORM does so when the design point lies beyond them.

PySTRA evaluates each tail on the side where its probability is small. With the survival function \(\bar F_X = 1 - F_X\),

\[\begin{split}x = \begin{cases} F_X^{-1}\bigl(\Phi(u)\bigr), & u \le 3,\\ \bar F_X^{-1}\bigl(\Phi(-u)\bigr), & u > 3, \end{cases} \qquad u = \begin{cases} \Phi^{-1}\bigl(F_X(x)\bigr), & F_X(x) \le \Phi(3),\\ -\Phi^{-1}\bigl(\bar F_X(x)\bigr), & \text{otherwise}. \end{cases}\end{split}\]

Up to \(u = 3\) the complement \(1 - \Phi(u)\) keeps its relative precision to about \(10^{-13}\). Both forms stay exact while the tail probability is a normal double, down to about \(2.2 \times 10^{-308}\) at \(|u| \approx 37.52\). PySTRA switches there to inverting the logarithm \(\log \Phi(u)\), before subnormal probabilities lose precision; they remain representable to \(|u| \approx 38.5\), although a CDF implementation may underflow earlier. The logarithm is inverted in closed form for the normal, lognormal and Gumbel families, and otherwise by solving \(\log F_X(x) = \log \Phi(u)\). Quantities such as \(\log(1 - e^{a})\) are computed as \(\log(-\operatorname{expm1} a)\) or \(\operatorname{log1p}(-e^{a})\), whichever avoids cancellation [Maechler2012].

Composite distributions inherit this accuracy from their parents. The maximum of \(N\) draws has \(\log F = N \log F_P\), and its survival probability \(q\) corresponds to the parent survival probability \(-\operatorname{expm1}\bigl(\operatorname{log1p}(-q)/N\bigr)\), which equals \(q/N\) to double precision once \(q < 10^{-200}\).

Numerical inverse CDFs can fail silently far into a tail: SciPy’s beta quantile stalls near \(4 \times 10^{-50}\) for probabilities below about \(10^{-100}\). PySTRA checks every quantile whose tail probability is below \(10^{-8}\) against the log-CDF and solves again where it misses. Jacobians \(f_X(x)/\varphi(u)\) are formed from log densities where either density underflows.

Two limits remain. A bounded support cannot place \(x\) closer to its bound than one unit in the last place, so a bounded tail loses resolution before its probability does. The conditional distributions of non-elliptical copulas, such as Frank, are still evaluated in probability space. Working at high reliability shows these properties in use.

Use this method: Copulas and isoprobabilistic transformations · Dependence and probability transformations · Distributions and dependence

For coordinate conventions, see Notation and probability spaces.