monte_carlo_double_ended draws non-section alpha with the wrong joint distribution

[bdestombe]

Two distinct silent bugs in the parameter draw block of monte_carlo_double_ended (src/dtscalibration/dts_accessor.py:1748-1766). Both attack the same code path. Issue #230 (RNG reproducibility) touches the same area; this issue is about the joint distribution of the draws, not the seed.

1. np.any(not_ix_sec) returns False when the only non-section index is 0

Location: src/dtscalibration/dts_accessor.py:1754

not_ix_sec = np.array([i for i in range(no) if i not in ix_sec])

if np.any(not_ix_sec):     # bug
    ...
    alpha[:, not_ix_sec] = not_alpha_mc

np.any tests element truthiness. np.any(np.array([0])) is False because 0 is falsy. When the only non-section position is x = 0 (e.g. the first reference section starts at the second pixel and covers everything beyond), the entire alpha-MC block for non-section x is skipped. alpha[:, 0] stays at zero across all MC realizations; tmpf_mc_set[:, 0, :] and tmpb_mc_set[:, 0, :] use a fixed-zero alpha rather than draws. CIs at x=0 understate the alpha-uncertainty contribution.

Fix: if not_ix_sec.size > 0:.

2. Non-section alpha is sampled from independent univariate Gaussians using diagonal-only variance

Location: src/dtscalibration/dts_accessor.py:1758-1762

not_alpha_var = p_cov[2*nt + 1 + not_ix_sec, 2*nt + 1 + not_ix_sec]   # diagonal only

not_alpha_mc = np.random.normal(
    loc=not_alpha_val,
    scale=not_alpha_var**0.5,
    size=(mc_sample_size, not_alpha_val.size),
)                                                                       # univariate, independent per x

Section-alpha is sampled jointly with (γ, df, db) via multivariate_normal.rvs(po_cov) (line 1736). Non-section-alpha is sampled outside that joint draw with only its diagonal variance. Result: in the MC sample,

• Cov(α_x_nonsec, γ) is 0 (true value generically nonzero),
• Cov(α_x_nonsec, α_x_sec) is 0 (true value generically nonzero),
• Cov(α_x_nonsec, α_x'_nonsec) is 0 for distinct non-section x (true value generically nonzero).

Per-(x,t) tmpf_mc_var is unaffected (correct mean and marginal variance). Aggregate quantities — spatial means, IVW pools, _avgx2, x-marginalized percentiles — are biased: they are computed under an independence assumption that the underlying joint posterior does not satisfy.

Fix: extend from_i (line 1725) to include all alpha indices np.arange(1+2*nt, 1+2*nt+no) so po_mc directly samples the full joint via multivariate_normal. The entire not_ix_sec block (and its np.any guard) can then be deleted. If memory profile shows the full-rank draw is too large, the alternative is to draw non-section alpha conditional on the section sample, but that is substantially more code.

Reproducer / tests

• Build a dataset with sections that omit only x = 0 (so not_ix_sec.tolist() == [0]); run monte_carlo_double_ended; assert the variance of alpha_mc[:, 0] across MC realizations is non-zero. Currently zero.
• Build a dataset with np.any(not_ix_sec) == True and check np.std(alpha_mc[:, x]) against sqrt(p_cov.diag()[2*nt+1+x]) (correct), then check np.cov(alpha_mc[:, [x1, x2]]) against the corresponding p_cov block (currently diag, should be full block).