Fix output-axis reductions and DR variance propagation

[kbattocchi]

Summary

Following up on #1067, had Copilot audit the codebase for other places where we wrongly use a numpy operation on a whole array instead of specifying an axis.

This PR fixes the four confirmed bugs that the review found, in three logical commits:

• Normalize output-averaged CV MSE by mean within-output target variance, rather than pooling output means into the variance. Independent outcome offsets no longer spuriously inflate model-selection scores.
• Scope population Gaussian-mixture interval invalidation to individual outcome/treatment columns. Exclude nonfinite means/scales and nonpositive scales before solving, preserving unaffected intervals and avoiding nontermination.
• Align DR sample_var vectors and singleton columns with final targets, avoiding unintended sample-by-sample broadcasting and fitting failures.
• Propagate DR outcome variance separately for each treatment-versus-control contrast using its squared observed-outcome derivative. Observations from unrelated treatment arms contribute zero propagated variance rather than inflating standard errors.

Includes regression tests and directly related documentation. Both separate and multitask DR final models are covered, including trimming, heterogeneous weights/propensities, and summarized-versus-expanded observation equivalence.