Copulas and isoprobabilistic transformations#

For worked examples, see Dependence and probability transformations; for the derivations, see Copulas as joint distribution specifications.

A joint distribution has two separate ingredients: marginal distributions and a copula describing their dependence. Choosing a transformation for reliability analysis does not change that joint distribution. This follows the distinction developed by Lebrun and Dutfoy [LebrunDutfoy2009a] [LebrunDutfoy2009b].

Specifying a joint distribution#

import pystra as ra

joint = ra.JointDistribution(
    [ra.Lognormal("R", 10.0, 2.0), ra.Normal("S", 4.0, 1.0)],
    ra.StudentTCopula([[1.0, 0.4], [0.4, 1.0]], df=4),
)
model = ra.StochasticModel(joint)
model.add_variable(ra.Constant("C", 1.0))
samples = joint.rvs(1000, seed=2026)  # rows of observations, columns R, S
print(joint.pdf([10.0, 4.0]))
print(joint.cdf([10.0, 4.0]))

Alternatively, add all random variables to an existing model and then call model.set_copula(copula). Constants may be added afterwards. The copula’s dimension must equal the number of random variables. Named continuous PySTRA marginals, including ScipyDist, are supported; mixed/discrete marginals are outside this implementation’s scope. Known zero-inflated marginals with positive point mass are rejected.

The initial families are:

  • IndependentCopula(dimension): the product copula.

  • GaussianCopula(R): R is the latent normal correlation.

  • StudentTCopula(R, df): R is the latent Student-t shape matrix with unit diagonal. Positive degrees of freedom are supported, including values for which the latent representative has no covariance. Identity R does not give independence, because all coordinates share a random scale.

  • FrankCopula(theta): a bivariate non-elliptical example; zero gives independence. The initial numerical implementation supports abs(theta) <= 30.

Gaussian and Student-t matrices must be symmetric and positive definite. They can also be constructed using from_kendall_tau(tau, ...), with df supplied for Student-t. This uses \(R_{ij}=\sin(\pi\tau_{ij}/2)\) and validates the resulting matrix. These matrices generally differ from the physical marginals’ Pearson correlations. Copula fitting and physical Pearson calibration for Frank or Student-t copulas are not included.

Existing model.set_correlation(...) retains its physical Pearson meaning and the existing Gaussian Nataf calibration. Setting a copula explicitly replaces that specification; setting a correlation explicitly switches back to the legacy specification. get_correlation() raises for an explicit copula to prevent confusion between physical and latent correlations. Use get_copula() for the dependence specification. get_joint_distribution() also works for a legacy model by calibrating its Gaussian copula.

Choosing the transformation#

# Condition on S before R
options = ra.FORMOptions(transform="rosenblatt", rosenblatt_order=[1, 0])
form = ra.FORM(
    model, ra.LimitState(lambda R, S, C: R - C*S), options=options
)
form_result = form.run()
print(form_result.design_index, form_result.failure_probability)

Rosenblatt uses sequential conditional CDFs followed by the normal quantile function. Its standard coordinates are independent standard normals for every supported copula. Returned coordinate arrays retain original variable indices: u[order[k]] is the k-th conditional innovation. Every component of a system must use the same order.

With default analysis options, an explicit Gaussian copula uses Nataf; Student-t and Frank use Rosenblatt. Existing models without an explicit copula retain their previous Gaussian Nataf behavior.

FORMOptions(transform="nataf") selects generalized Nataf for an elliptical copula. "cholesky" and "svd" select Nataf factorizations explicitly; these choices are unavailable for Frank. Conditioning order applies only to Rosenblatt. Gaussian Nataf with Cholesky factorization equals Rosenblatt in the same order. Different Gaussian orders are related by an orthogonal change of normal coordinates, so the optimum FORM probability is invariant. For non-Gaussian copulas, changing Rosenblatt order can change the FORM approximation, although it cannot change the exact failure probability.

Transformations are also available independently of an analysis:

transform = joint.make_transformation("rosenblatt", order=[1, 0])
u = transform.x_to_u([10.0, 4.0])
x = transform.u_to_x(u)
du_dx = transform.jacobian_u_wrt_x(u, x)

These methods take one full vector; joint distribution methods take either one point or rows of points. Copula-level rosenblatt and inverse_rosenblatt operate on uniform coordinates, whereas the joint transformation uses normal coordinates. Elliptical transformation Jacobians are analytic; Frank uses central differences of the inverse transformation, without extra limit-state evaluations.

Generalized Nataf is not always normal#

For an elliptical copula, generalized Nataf first maps each marginal into the elliptical representative and then removes its shape matrix. The result is spherical, but is independent normal only for a Gaussian copula.

For StudentTCopula(R, df), the Nataf coordinates have spherical Student-t law with identity shape, not independent coordinates or unit covariance. The univariate scale is one; when df > 2, covariance is \(\nu/(\nu-2) I\). FORM therefore uses

\[P_{f,\mathrm{FORM}} = T_\nu(-\beta),\]

where \(\beta\) is the signed geometric distance in this spherical space. form_result.design_index is this distance, while form_result.beta is \(-\Phi^{-1}(P_f)\) for comparison with conventional normal reliability indices.

SORM, system FORM, simulation methods and the Strong Maximum Test currently require independent normal coordinates and reject spherical Student-t space. Use Rosenblatt with these methods. Closed-form parameter sensitivities still assume the legacy physical-Pearson Nataf specification; for explicit copulas, use numerical sensitivities, which keep the copula parameters fixed.

Validation and numerical scope#

The tests compare Gaussian and Student-t joint densities with SciPy’s multivariate distributions, verify transformation Jacobians against density identities, transform independent reference Student-t samples, and check an exact affine Student-t FORM half-space probability.

The second paper’s exponential benchmark uses rates 1 and 3 and failure event \(8X_1+2X_2-1\leq0\). With Frank parameter 10, the two Rosenblatt orders reproduce FORM probabilities approximately 0.107 and 0.122; direct integration gives approximately 0.1038. With Gaussian latent correlation 0.5, both orders and Cholesky Nataf give approximately 0.09758.

System reliability and conditioning order#

Meinen and Steenbergen (2025), Application of the Rosenblatt transformation in First-Order System Reliability approximations, extend the exponential example to a two-component series system. Example 1 uses Gaussian latent correlation 0.5 or Frank parameter 4.73, rather than the earlier benchmark’s 10. For fully dependent identical components, both components use the same pair of variables, so their union must equal a single component event.

The worked tutorial credits Meinen and Steenbergen and develops their example with transformation plots, shared-event and independent-copy systems, coordinate alignment, integration and direct simulation. tests/test_rosenblatt_system_order.py preserves the shared-order results as regressions:

Fully dependent identical components: system failure probability#

Copula

Canonical Rosenblatt / FORM

Reverse Rosenblatt / FORM

Direct original-event integration

Gaussian, latent correlation 0.5

0.0975813

0.0975813

0.0871846

Frank, theta 4.73

0.1019922

0.1142769

0.0932253

The direct reference integrates the physical failure event along each coordinate separately; both evaluations agree. The dependent-system CMC probabilities in Section 3.3 of the paper appear transposed between Gaussian and Frank. The values above use the stated copula parameters and marginal rates, supported by direct sampling, rather than adopting that apparent transcription error.

A fully specified copula fixes the joint law. It does not make the FORM approximation invariant to a nonlinear change of reference coordinates. Gaussian Rosenblatt orders are related by an orthogonal transformation: applying one common order to every component preserves the tangent-event correlations and system probability. For a Frank copula, changing order can change the transformed failure-surface geometry and hence the approximation. The original-event probability remains independent of either choice.

Alpha vectors from different conditioning orders cannot simply be dotted together, even when their entries carry the same original variable names. Deliberately doing so reproduces the paper’s invalid mixed-order probabilities of approximately 0.133 (Gaussian) and 0.156 (Frank) for identical events. SystemFORM avoids this construction by analyzing every component against one complete model and one transformation configuration. A second regression uses the paper’s distinct-component variation, verifying Gaussian system invariance between Nataf and both Rosenblatt orders.

These are first-order system approximations. Their normal-score dependence is still represented by the tangent-plane alpha products; a non-Gaussian input copula does not turn those planes into an exact system failure event. Use original-event integration or simulation to assess approximation error. Exploring several orders is a sensitivity diagnostic, not a certified bound or an automatic prescription to select the largest probability.

The same three regressions pass against the v1 copula feature commit f3d0d1f after adapting the four changed API spellings. They are suitable for backport with that feature; this does not imply copula support is present in the released v1 package. The paper’s larger Example 2 also needs explicit joint dependence across components. Its local-component autocorrelation approximation should not be substituted silently for a complete joint law. PySTRA’s Frank implementation currently supports two variables; composition of independent Frank blocks is a separate extension.

Numerical limits#

Transformations require finite points and resolvable interior probabilities; they do not silently clip a failed tail calculation into the unit interval. Gaussian and Student-t joint transformations use survival functions where available to preserve tail resolution. Copula uniform coordinates still have floating-point limits: in very low conditional-density regions, roundoff can be amplified substantially by the inverse map. Marginal implementations may impose additional tail limits. Joint CDFs for elliptical copulas use numerical SciPy integration and accept its integration keyword arguments; they are not exact rare-event probability estimators.

Continue: Dependence and probability transformations · Distributions and dependence · Probability transformations and copulas