1import numpy as np2from scipy import stats34from chronos_x.eval.scores import crps_student_t, crps_gaussian, crps_samples, weighted_interval_score, coverage, energy_score5from chronos_x.eval.statistics import diebold_mariano, cluster_bootstrap_mean, holm, benjamini_hochberg, hac_variance678def test_crps_student_t_matches_monte_carlo():9 rng = np.random.default_rng(0)10 y = np.array([-0.001, 0.0, 0.002, 0.004]); mu = np.array([0.0002, 0.0, -0.0001, 0.0]); sc = np.array([0.001, 0.002, 0.0008, 0.0015]); nu = np.array([3.0, 4.5, 8.0, 30.0])11 cf = crps_student_t(y, mu, sc, nu)12 s = stats.t.rvs(nu[None, :], loc=mu[None, :], scale=sc[None, :], size=(200_000, 4), random_state=rng)13 mc = crps_samples(y, s)14 np.testing.assert_allclose(cf, mc, rtol=0.02)151617def test_crps_student_t_large_nu_tends_to_gaussian():18 y = np.array([0.001]); mu = np.array([0.0]); sc = np.array([0.0012])19 assert abs(crps_student_t(y, mu, sc, np.array([1e6]))[0] - crps_gaussian(y, mu, sc)[0]) < 1e-9202122def test_crps_samples_degenerate_equals_abs_error():23 y = np.array([0.5, -1.0]); s = np.tile(np.array([0.0, 0.0]), (50, 1))24 np.testing.assert_allclose(crps_samples(y, s), np.abs(y), atol=1e-12)252627def test_wis_and_coverage_sane():28 rng = np.random.default_rng(1)29 s = rng.normal(0, 1, (2000, 500)); y = rng.normal(0, 1, 500)30 c, w = coverage(y, s, 0.9)31 assert 0.85 < c < 0.95 and 3.0 < w < 3.632 assert 0.3 < weighted_interval_score(y, s).mean() < 0.8333435def test_energy_score_prefers_true_dependence():36 rng = np.random.default_rng(2)37 N, H, n = 400, 10, 30038 L = np.linalg.cholesky(0.6 * np.ones((H, H)) + 0.4 * np.eye(H))39 Y = rng.normal(size=(N, H)) @ L.T40 S_dep = rng.normal(size=(n, N, H)) @ L.T41 S_ind = rng.normal(size=(n, N, H))42 assert energy_score(Y, S_dep).mean() < energy_score(Y, S_ind).mean()434445def test_dm_hac_and_bootstrap_and_multiplicity():46 rng = np.random.default_rng(3)47 n = 300048 e = rng.normal(size=n); d = e + 0.049 # MA(4)-overlap style autocorrelation50 d = np.convolve(e, np.ones(5) / 5, mode="same")51 r = diebold_mariano(d, np.zeros(n), np.arange(n), lag=4)52 assert abs(r["t"]) < 453 assert hac_variance(d, 4) > d.var() * 2 # positive autocorrelation inflates the long-run variance54 cb = cluster_bootstrap_mean(d, np.repeat(np.arange(n // 50), 50), n_boot=300)55 assert cb["lo"] <= cb["mean"] <= cb["hi"]56 assert np.allclose(holm([0.01, 0.04, 0.03]), [0.03, 0.06, 0.06])57 assert np.allclose(benjamini_hochberg([0.01, 0.04, 0.03]), [0.03, 0.04, 0.04])58