Binary / non-Gaussian-outcome IV via control-function compilation

[daimon-pymclabs]

Part of #343. Depends on the IV-aware identification issue; builds on the continuous-outcome example.

Problem

The literal binary-outcome IV case (e.g. a Bernoulli/logit purchase outcome with an endogenous continuous price) cannot currently be expressed via ~~: _validate_residual_cov_families in compile.py hard-requires every variable in a residual-covariance block to be Gaussian, because the block compiles to a joint MvNormal. A Bernoulli outcome has no Gaussian residual to correlate.

Proposal

Add a control-function compilation path for residual-covariance blocks that mix a non-Gaussian outcome with a Gaussian treatment equation:

• keep the treatment equation Gaussian (treatment ~ instrument + controls);
• expose its residual as a Deterministic;
• inject that residual into the non-Gaussian outcome's linear predictor (e.g. the logit) with a learned coefficient.

The compiler selects control-function vs joint-MvNormal based on the families present in the block — so this needs no user-facing API change: the user still writes treatment ~~ outcome + the instrument and calls fit(). Relax / branch _validate_residual_cov_families accordingly.

Why this is the Bayesian win

Because the first-stage residual is a deterministic function of the sampled stage-1 parameters within one joint model, second-stage uncertainty propagates correctly for free — no bootstrap or analytic SE correction, unlike the frequentist two-step control function.

Acceptance

• A gaussian-treatment + bernoulli-outcome block compiles via the control-function path.
• Recovers a known causal coefficient on simulated data with correct posterior coverage.
• Existing Gaussian-only ~~ behaviour unchanged.