1import json
  2from pathlib import Path
  3import numpy as np
  4from scipy import stats, optimize
  5
  6RNG = np.random.default_rng(1434)
  7EPS = 1e-10
  8
  9
 10def copula_logpdf(u1, u2, rho):
 11    """log h1(u1) + log h2(u2|u1) for a Gaussian copula."""
 12    z1, z2 = stats.norm.ppf(np.clip(u1, EPS, 1-EPS)), stats.norm.ppf(np.clip(u2, EPS, 1-EPS))
 13    rho = float(np.clip(rho, -0.999999, 0.999999))
 14    return (-0.5*np.log1p(-rho*rho)
 15            - (rho*rho*(z1*z1 + z2*z2) - 2*rho*z1*z2)/(2*(1-rho*rho)))
 16
 17
 18def conditional_cdf(u2, u1, rho):
 19    z1, z2 = stats.norm.ppf(np.clip(u1, EPS, 1-EPS)), stats.norm.ppf(np.clip(u2, EPS, 1-EPS))
 20    return stats.norm.cdf((z2-rho*z1)/np.sqrt(1-rho*rho))
 21
 22
 23def make_sample(n, rho, marginal_pair=("lognormal", "gamma"), seed=0):
 24    rng = np.random.default_rng(seed)
 25    z = rng.normal(size=(n, 2))
 26    z[:, 1] = rho*z[:, 0] + np.sqrt(1-rho*rho)*z[:, 1]
 27    u = stats.norm.cdf(z)
 28    if marginal_pair[0] == "lognormal":
 29        x1 = stats.lognorm.ppf(u[:, 0], s=0.8, scale=np.exp(0.3))
 30        f1 = stats.lognorm.pdf(x1, s=0.8, scale=np.exp(0.3))
 31    else:
 32        x1 = stats.gamma.ppf(u[:, 0], a=2.0, scale=1.2)
 33        f1 = stats.gamma.pdf(x1, a=2.0, scale=1.2)
 34    if marginal_pair[1] == "gamma":
 35        x2 = stats.gamma.ppf(u[:, 1], a=3.0, scale=0.7)
 36        f2 = stats.gamma.pdf(x2, a=3.0, scale=0.7)
 37    else:
 38        x2 = stats.lognorm.ppf(u[:, 1], s=0.5, scale=np.exp(0.2))
 39        f2 = stats.lognorm.pdf(x2, s=0.5, scale=np.exp(0.2))
 40    return x1, x2, u, f1, f2
 41
 42
 43def fit_rho(u):
 44    def objective(a):
 45        return -np.mean(copula_logpdf(u[:, 0], u[:, 1], np.tanh(a)))
 46    result = optimize.minimize_scalar(objective, bounds=(-4, 4), method="bounded")
 47    return float(np.tanh(result.x))
 48
 49
 50def pit_sweep():
 51    rows = []
 52    settings = [("lognormal", "gamma"), ("gamma", "lognormal"), ("lognormal", "lognormal"), ("gamma", "gamma")]
 53    for j, pair in enumerate(settings):
 54        _, _, u, _, _ = make_sample(30000, 0.0, pair, 100+j)
 55        for k in range(2):
 56            rows.append({"marginals": "%s/%s" % pair, "coordinate": k+1,
 57                         "mean": float(np.mean(u[:, k])), "variance": float(np.var(u[:, k])),
 58                         "ks": float(stats.kstest(u[:, k], "uniform").statistic)})
 59    return rows
 60
 61
 62def dependence_sweep():
 63    rows = []
 64    for i, rho in enumerate([0.0, 0.2, 0.5, 0.8]):
 65        _, _, u, _, _ = make_sample(40000, rho, ("lognormal", "gamma"), 200+i)
 66        observed = float(stats.spearmanr(u[:, 0], u[:, 1]).statistic)
 67        predicted = float(6/np.pi*np.arcsin(rho/2))
 68        kendall_observed = float(stats.kendalltau(u[:, 0], u[:, 1]).statistic)
 69        kendall_predicted = float(2/np.pi*np.arcsin(rho))
 70        rows.append({"rho": rho, "predicted_spearman": predicted, "observed_spearman": observed,
 71                     "abs_error": abs(observed-predicted),
 72                     "predicted_kendall": kendall_predicted,
 73                     "observed_kendall": kendall_observed,
 74                     "kendall_abs_error": abs(kendall_observed-kendall_predicted)})
 75    return rows
 76
 77
 78def likelihood_sweep():
 79    rows = []
 80    for i, rho in enumerate([0.0, 0.2, 0.5, 0.8]):
 81        _, _, u, f1, f2 = make_sample(50000, rho, ("lognormal", "gamma"), 300+i)
 82        marginal_nll = float(np.mean(-np.log(np.maximum(f1, EPS))-np.log(np.maximum(f2, EPS))))
 83        indep = marginal_nll
 84        joint = float(np.mean(-np.log(np.maximum(f1, EPS))-np.log(np.maximum(f2, EPS))
 85                              -copula_logpdf(u[:, 0], u[:, 1], rho)))
 86        observed_gain = indep-joint
 87        predicted_gain = float(-0.5*np.log(1-rho*rho))
 88        fitted = fit_rho(u)
 89        rows.append({"rho": rho, "independent_nll": indep, "copula_nll": joint,
 90                     "observed_nll_gain": observed_gain, "predicted_nll_gain": predicted_gain,
 91                     "fitted_rho": fitted})
 92    return rows
 93
 94
 95def derivative_check():
 96    u1, u2, rho = 0.37, 0.61, 0.65
 97    h = 1e-5
 98    numerical = (conditional_cdf(u2+h, u1, rho)-conditional_cdf(u2-h, u1, rho))/(2*h)
 99    z1, z2 = stats.norm.ppf(u1), stats.norm.ppf(u2)
100    analytic = np.exp(copula_logpdf(np.array([u1]), np.array([u2]), rho))[0]
101    return {"u1": u1, "u2": u2, "rho": rho, "finite_difference": float(numerical),
102            "analytic_density": float(analytic), "relative_error": float(abs(numerical-analytic)/analytic)}
103
104
105def main():
106    pit = pit_sweep()
107    dep = dependence_sweep()
108    lik = likelihood_sweep()
109    deriv = derivative_check()
110    # A direct baseline-vs-idea mini comparison at substantial dependence.
111    target = lik[-1]
112    summary = {
113        "predictions": {
114            "PIT_uniformity": "mean -> 0.5, variance -> 1/12, and KS remains small for every continuous marginal",
115            "rank_invariance": "For a Gaussian copula, Spearman = (6/pi) asin(rho/2) and Kendall tau = (2/pi) asin(rho), independent of marginal family",
116            "dependence_gain": "independent NLL - copula NLL = -0.5 log(1-rho^2), exactly zero at rho=0"
117        },
118        "derivative_check": deriv,
119        "pit_sweep": pit,
120        "dependence_sweep": dep,
121        "likelihood_sweep": lik,
122        "mini_comparison": {"dependence": target["rho"], "baseline_independent_nll": target["independent_nll"],
123                            "idea_copula_nll": target["copula_nll"], "idea_fitted_rho": target["fitted_rho"]}
124    }
125    Path("results.json").write_text(json.dumps(summary, indent=2))
126    print(json.dumps(summary, indent=2))
127
128if __name__ == "__main__":
129    main()