Preserving Copula Correlations in Synthetic Attributes

When you generate each attribute of a synthetic parcel or sensor record from its own marginal distribution, the joint dependence — population correlating with built area, elevation anti-correlating with flood risk — vanishes, and a Gaussian copula is the mechanism that restores it without distorting any single attribute’s shape.

Part of Attribute Correlation Modeling for Synthetic Spatial Data: that page covers dependence modeling across attribute types in general; this page drills into the specific, dominant technique — the Gaussian copula — showing how to fit marginals independently, inject the rank-correlation structure through a Gaussian dependence layer, and gate the result on a correlation-matrix distance so a decorrelated batch never reaches a consumer.

Root Cause: Why Independent Marginals Destroy Joint Structure

A copula separates a joint distribution into two parts: the shapes of the individual variables (their marginals) and the way they move together (their dependence structure). Sklar’s theorem guarantees that any multivariate distribution FF can be written as

F(x1,…,xd)=C(F1(x1),…,Fd(xd))F(x_1, \dots, x_d) = C\big(F_1(x_1), \dots, F_d(x_d)\big)

where each FiF_i is a marginal CDF and CC is a copula — a distribution on the unit hypercube [0,1]d[0,1]^d with uniform margins that carries all of the dependence. The practical consequence is stark: if you sample each attribute independently from its fitted marginal, you have implicitly chosen the independence copula C(u1,…,ud)=∏uiC(u_1,\dots,u_d)=\prod u_i, which throws away every pairwise and higher-order relationship. The marginals can each be pixel-perfect and the joint distribution still be wrong.

For spatial attributes this failure is quietly destructive because downstream models lean on the correlations, not the marginals. A flood-risk model trained on synthetic parcels where elevation and flood depth are independent learns nothing useful; a demand model trained where population and road density are decorrelated produces confident nonsense. The marginals pass every univariate check — histograms match, means and variances match — so the corruption survives naive validation and only surfaces as degraded model accuracy downstream.

The Gaussian copula reproduces dependence through rank correlation rather than linear (Pearson) correlation, which matters for spatial attributes that are heavy-tailed or skewed. Rank correlation (Spearman’s ρ\rho, Kendall’s τ\tau) is invariant to any monotone transform of a marginal, so it survives the log-transforms and power-transforms that skewed spatial variables routinely need. You fit the dependence in a rank/normal-score space where it is stable, then map back through each marginal’s inverse CDF. The correlations that must be preserved are the same ones the realism metrics evaluation stage scores, so the copula and the gate share one definition of correctness.

Gaussian copula synthesis of correlated spatial attributes with a correlation-matrix-distance gate Top row, the fit path left to right: real attributes X (population, built area, elevation) pass through their empirical CDFs to uniform ranks U, then through the inverse standard-normal to normal scores Z; a correlation matrix R is estimated from Z. Bottom row, the sample path right to left: correlated Gaussian samples are drawn from N(0, R), mapped through the standard-normal CDF back to uniforms U-tilde, then through independently fitted inverse marginals to produce synthetic attributes X-tilde whose marginals match the originals but now carry the dependence. A gate box compares the correlation matrix of X-tilde against R using a matrix-distance metric and blocks the batch if the distance exceeds a tolerance. Fit path → estimate dependence Real X pop, built, elev Ranks U empirical CDF Scores Z Φ⁻¹(U) Correlation R from normal scores Sample path ← reproduce dependence Draw N(0, R) correlated Gaussians Uniforms Ũ Φ(Z̃) Inverse Fᵢ⁻¹ fitted per attribute Synthetic X̃ marginals + dependence Correlation-matrix-distance gate d(corr(X̃), R) ≤ tol → promote · else block
Dependence is estimated once as a correlation matrix R in normal-score space, then replayed into synthetic attributes whose marginals are fitted independently; the gate blocks any batch whose correlation matrix drifts from R.

Minimal Reproducer: Show Independent Sampling Kills Correlation

Pin the toolchain, then sample each attribute from its correct marginal independently and watch the real correlation collapse toward zero.

python
# requirements.txt — pin majors; let patch releases float
numpy==1.26.*
scipy==1.11.*

python
import numpy as np
from scipy import stats

rng = np.random.default_rng(11)
# Real attributes with strong positive dependence (population vs built area).
z = rng.multivariate_normal([0, 0], [[1, 0.82], [0.82, 1]], size=5000)
pop = np.expm1(z[:, 0] * 0.6 + 3)          # skewed, log-normal-ish marginal
built = np.expm1(z[:, 1] * 0.5 + 2)
real_rho = stats.spearmanr(pop, built).statistic          # ~0.80

# Naive "synthetic": resample each marginal INDEPENDENTLY -> dependence destroyed.
syn_pop = rng.permutation(pop)
syn_built = rng.permutation(built)
naive_rho = stats.spearmanr(syn_pop, syn_built).statistic  # ~0.00
print(round(real_rho, 2), round(naive_rho, 2))             # 0.80 0.0

The marginals of syn_pop and syn_built are identical to the real ones — every univariate check passes — yet the rank correlation has gone to zero. That is the failure the copula fixes.

Fix: Fit Marginals Independently, Inject a Gaussian Dependence Layer

The construction has three parts: estimate the copula correlation from normal scores, fit each marginal separately, then sample correlated Gaussians and push them through the inverse marginals.

Step 1 — Estimate the copula correlation in normal-score space

Convert each attribute to uniform ranks with its empirical CDF, then to standard-normal scores, and estimate the correlation there. Working in rank space makes the estimate invariant to the skew of the raw marginals.

python
import numpy as np
from scipy import stats

def fit_copula_correlation(X: np.ndarray) -> np.ndarray:
    """Correlation matrix R of the Gaussian copula, estimated from normal scores.
    X: (n_samples, n_attributes) real attribute matrix."""
    n, d = X.shape
    ranks = np.empty_like(X, dtype="float64")
    for j in range(d):
        # average ranks -> uniforms in (0, 1), open interval avoids +/-inf below
        ranks[:, j] = stats.rankdata(X[:, j], method="average") / (n + 1)
    z = stats.norm.ppf(ranks)                 # normal scores
    R = np.corrcoef(z, rowvar=False)
    return R

Step 2 — Fit each marginal independently

Fit or store each attribute’s marginal on its own. Independent marginal fitting is a feature, not a shortcut: it lets you model a heavy-tailed population count and a bounded slope percentage with entirely different families while the copula handles their joint behaviour. Use the empirical inverse CDF (a quantile function) to avoid forcing a parametric family that does not fit.

python
import numpy as np

class EmpiricalMarginal:
    """Invertible empirical marginal: uniform u in (0,1) -> attribute value."""
    def __init__(self, x: np.ndarray):
        self.sorted = np.sort(x.astype("float64"))
    def ppf(self, u: np.ndarray) -> np.ndarray:
        # linear-interpolated empirical quantile function
        idx = u * (len(self.sorted) - 1)
        lo = np.floor(idx).astype(int)
        frac = idx - lo
        hi = np.minimum(lo + 1, len(self.sorted) - 1)
        return self.sorted[lo] * (1 - frac) + self.sorted[hi] * frac

def fit_marginals(X: np.ndarray) -> list[EmpiricalMarginal]:
    return [EmpiricalMarginal(X[:, j]) for j in range(X.shape[1])]

Step 3 — Sample correlated attributes from a seeded stream

Draw correlated standard normals via the Cholesky factor of RR, map to uniforms through the normal CDF, then through each inverse marginal. Seed from an explicit integer so the synthetic batch is reproducible, matching the reproducibility contract the point process simulation models hold generated geometry to.

python
import numpy as np
from scipy import stats

def sample_correlated(R: np.ndarray, marginals: list[EmpiricalMarginal],
                      n: int, *, seed: int) -> np.ndarray:
    """Draw n synthetic rows carrying the copula dependence R and the given marginals."""
    rng = np.random.default_rng(seed)
    d = R.shape[0]
    # Repair R to the nearest valid correlation matrix if needed (see gotchas).
    L = np.linalg.cholesky(R)
    z = rng.standard_normal((n, d)) @ L.T          # correlated normals, corr -> R
    u = stats.norm.cdf(z)                          # uniforms carrying the dependence
    return np.column_stack([m.ppf(u[:, j]) for j, m in enumerate(marginals)])

Because the dependence lives entirely in RR and the shapes live entirely in the marginals, you can swap either independently — regenerate with a stronger correlation without refitting marginals, or update one marginal without touching the dependence. This clean separation is why the copula is the default over jointly training a GAN-based spatial generation model when only pairwise dependence matters.

Verification Step: Gate on a Correlation-Matrix Distance

Gate the synthetic batch on how far its correlation matrix has drifted from the fitted RR. A single scalar distance over the whole matrix catches decorrelation that per-attribute marginal checks miss entirely.

python
import numpy as np
from scipy import stats

def corr_matrix_distance(A: np.ndarray, B: np.ndarray) -> float:
    """Correlation-matrix distance in [0, 1]: 0 == identical structure."""
    fa, fb = A.ravel(), B.ravel()
    return float(1.0 - (fa @ fb) / (np.linalg.norm(fa) * np.linalg.norm(fb)))

def test_copula_preserves_dependence(X_real, X_syn, R, tol=0.02):
    # 1. Marginals must still match per attribute (KS on each column).
    for j in range(X_real.shape[1]):
        ks = stats.ks_2samp(X_real[:, j], X_syn[:, j]).statistic
        assert ks < 0.05, f"attribute {j} marginal drifted: KS={ks:.3f}"

    # 2. Rank-correlation structure must match R within tolerance.
    z_syn = stats.norm.ppf(
        np.column_stack([stats.rankdata(X_syn[:, j], method="average") for j in range(X_syn.shape[1])])
        / (len(X_syn) + 1))
    R_syn = np.corrcoef(z_syn, rowvar=False)
    d = corr_matrix_distance(R, R_syn)
    assert d < tol, f"dependence structure drifted: matrix distance {d:.4f} > {tol}"

The two assertions are complementary: the KS gate catches a marginal that has drifted while the correlation-matrix-distance gate catches decorrelation with intact marginals — exactly the failure the naive reproducer above produces and no univariate test can see.

Edge Cases & Gotchas

A near-singular correlation matrix breaks the Cholesky factor. When two attributes are almost perfectly correlated (elevation and negative flood depth, say), the estimated RR can be numerically non-positive-definite and np.linalg.cholesky raises. Project RR onto the nearest positive-definite correlation matrix — clip negative eigenvalues to a small floor and rescale the diagonal back to one — before factoring, and log the correction magnitude so a genuinely degenerate pair is not silently masked.

Ties in a discrete spatial attribute distort the ranks. Categorical or heavily-quantized attributes (zoning class, integer building counts) produce large blocks of tied values, so rankdata with average ranks maps many records to the same normal score and understates the copula correlation. Break ties with a tiny seeded jitter before ranking, or model discrete attributes with a dedicated discrete copula rather than forcing them through the Gaussian construction.

Tail dependence the Gaussian copula cannot represent. The Gaussian copula has zero tail dependence: it cannot reproduce the tendency of extreme values to co-occur (simultaneous extreme rainfall and extreme runoff across adjacent cells). If joint extremes matter for the consumer, the correlation-matrix gate will pass while the tails are wrong — validate the joint tail explicitly and switch to a t- or Gumbel copula when it fails, because no amount of tuning fixes a structural limitation of the family.

Frequently Asked Questions

Why use rank correlation instead of Pearson correlation for the copula?
Rank correlation is invariant to any monotone transform of a marginal, so it survives the log and power transforms that skewed spatial attributes routinely need, whereas Pearson correlation changes under those transforms and conflates dependence with marginal shape. Estimating the copula in rank or normal-score space gives a dependence estimate that is stable regardless of how heavy-tailed the raw attributes are, which is why the Gaussian copula is fitted there rather than on raw values.
Can I fit each marginal with a different distribution family?
Yes, and you should. The whole point of the copula construction is that marginals and dependence are separated, so you can model a heavy-tailed population count, a bounded percentage, and a discrete class with three different approaches while a single correlation matrix handles their joint behaviour. Using the empirical inverse CDF per attribute avoids forcing a parametric family that does not fit the data.
My marginals all pass their checks but the downstream model degraded — why?
Almost certainly the joint dependence was destroyed while the marginals stayed correct, which is exactly what independent sampling does. Univariate checks such as histograms and KS tests cannot see decorrelation because each attribute in isolation is perfect. Add a correlation-matrix-distance gate that compares the synthetic correlation matrix against the fitted one; it catches the failure that per-attribute checks structurally cannot.
When should I abandon the Gaussian copula for another family?
When joint extremes matter, because the Gaussian copula has zero tail dependence and cannot reproduce extreme values co-occurring across attributes or cells. If a validation of the joint tail fails while the correlation-matrix gate passes, switch to a t-copula for symmetric tail dependence or a Gumbel copula for upper-tail dependence. This is a structural limitation of the Gaussian family, not something you can tune away.