// gdpar Path 1 -- Empirical Bayes (EB) Marginal Template, Path B
// (Sub-phase 8.6.C, decision D34 canonized in Session 10 / 2026-05-25).
//
// Used by R/eb.R via cmdstanr::laplace() in Step (i) of the EB workflow
// for the regime K >= 2, p = 1 (multi-parametric distributional
// regression with univariate theta_ref per slot). Body is bit-identical
// to amm_distrib_K.stan; only the banner documents the EB context.
// Keeping the body bit-identical avoids subtle codegen drift between
// the FB and EB Path B paths and lets us reuse generate_stan_code_K()
// verbatim via the template_path override, including the full
// dispatch over stan_id in {1, 3, 5, 6, 7, 8, 9, 10, 11, 12, 13} that
// covers D33 deferral without additional work.
//
// The companion conditional template amm_eb_conditional_K.stan moves
// theta_ref_k from parameters{} to data{} (named theta_ref_k_data
// there), removes the per-slot hyperparameters mu_theta_ref_k /
// sigma_theta_ref_k from parameters{}, and removes the
// THETA_REF_PRIOR_BLOCK substitution from model{}, leaving the rest
// of the body bit-identical.
//
// Likelihood (Sub-decision 1 of Unit 3 paso 3, 2026-05-20; extended
// in sub-phase 8.3.4, 2026-05-20; extended in sub-phase 8.3.5a,
// 2026-05-21): distributional regression with canonical links per
// slot. For Gaussian K = 2:
//   y[i] ~ normal(eta_mu[i], exp(eta_sigma[i]))
// where eta_mu[i] follows the AMM canonical form on the identity-link
// scale with its own (theta_ref_mu, a_mu, b_mu) and
//   eta_sigma[i] = theta_ref_sigma + a_sigma(x_i)
//                  + b_sigma(x_i) * theta_ref_sigma
//                  + (W(theta_ref_sigma) - W(theta_anchor_sigma)) * x_i
// follows the AMM canonical form on the log-link scale (sigma_i =
// exp(eta_sigma[i])). For neg_binomial_2 K = 2 the analogue uses
// neg_binomial_2_log(eta_mu, exp(eta_phi)). For Beta K = 2 the
// canonical (mu, log phi) parametrization uses
// beta_proportion_lpdf(y_real | inv_logit(eta_mu), exp(eta_phi)).
// For Gamma K = 2 the canonical (mu, log shape) parametrization uses
// gamma_lpdf(y_real | exp(eta_shape), exp(eta_shape) / exp(eta_mu)).
// For Student-t K = 3 the canonical (mu, log sigma, log nu)
// parametrization uses
// student_t_lpdf(y_real | exp(eta_nu), eta_mu, exp(eta_sigma)).
// For Tweedie K = 3 the canonical (mu log, phi log, p identity)
// parametrization uses tweedie_lpdf(y_real | exp(eta_mu),
// exp(eta_phi), eta_p) with the hybrid Dunn--Smyth lpdf defined in
// the functions {} block.
// Built-in K-individual families wired in this template after
// sub-phase 8.3.6:
//   1  = gaussian              (mu identity, sigma log)              K=2
//   3  = neg_binomial_2        (mu log,      phi log)                K=2
//   5  = beta                  (mu logit,    phi log)                K=2
//   6  = gamma                 (mu log,      shape log)              K=2
//   7  = lognormal_loc_scale   (mu identity, sigma log) via D-A3.B   K=2
//   8  = student_t             (mu identity, sigma log,  nu log)     K=3
//   9  = tweedie               (mu log,      phi log,    p identity) K=3
//   10 = zip                   (mu log,      pi logit)               K=2
//   11 = zinb                  (mu log,      phi log,    pi logit)   K=3
//   12 = hurdle_poisson        (mu log,      pi logit)               K=2
//   13 = hurdle_neg_binomial_2 (mu log,      phi log,    pi logit)   K=3
// Continuous hurdle families (hurdle_gamma, hurdle_lognormal) and
// heterogeneous family/link per slot are deferred to sub-phases 8.3.7+.
// Stan-side custom-family wiring (D-A3.B) lives in 8.3.4 via the
// descriptor-driven registry .gdpar_K_custom_patterns().
//
// Scope of W (Sub-decision 2 of Unit 3 paso 3, 2026-05-20): W is
// globally shared across all K slots. A single W_raw, a single
// sigma_W and a single basis matrix feed the modulating term in
// every slot. Alternative scope (W per slot, or W active only in
// the location slot) is registered as a candidate decision for
// Session 8.4 if theoretical motivation arises.
//
// Scope of theta_anchor (Sub-decision 3 of Unit 3 paso 3,
// 2026-05-20): per-slot. Each K-individual slot carries its own
// theta_anchor_k[k]. Default 0. Stan data block exposes the vector
// theta_anchor_K of length K; the corresponding R-side input is a
// named numeric vector with names matching the slot names of the
// gdpar_formula_set (e.g. theta_anchor = c(mu = 0, sigma = 0)).
//
// Generated by R/stan_codegen.R via placeholder substitution; do
// not edit the substituted output by hand.
//
// Internal sampling parametrization. Same conventions as
// amm_distrib_multi.stan:
//   a    : configurable CP/NCP per slot via segment() priors. The
//          three a-block placeholders (DATA_CP_A_PER_K_DECL,
//          TP_A_BLOCK, MODEL_A_BLOCK) are resolved by
//          generate_a_blocks_K() from the cp_a_per_K logical vector
//          of length K.
//   W    : configurable CP/NCP globally (placeholders W_SCALE and
//          W_PRIOR).
//   b    : linear reparametrization per slot. Define c_b_k = theta_ref_k
//          * b_coef_k; the model samples c_b_k_raw directly (NCP on
//          c_b_k_raw) and reports b_coef_k as a derived quantity in
//          the single-anchor regime. This eliminates the bimodality
//          of (theta_ref_k, b_k) — see Recovery 2 (handoff 4 addendum,
//          2026-05-09).
//
// Centering of Z_a_k[k] and Z_b_k[k] is enforced by column-wise
// centering on the R side (see build_amm_design_K()); (C2)/(C3) of
// Block 1 hold identically per slot. (C4) holds per slot via
// W(theta_ref_k[k]) - W(theta_anchor_K[k]).

functions {
  // Tweedie compound Poisson-gamma log-density for the regime
  // 1 < p < 2 (sub-phase 8.3.5b of Block 8, 2026-05-21). Hybrid
  // dispatch per decision D2 of the scoping: Dunn--Smyth (2005)
  // infinite-series expansion in the central region
  // |p - 1.5| < tau (tau = 0.4 canonical) and the saddlepoint
  // approximation elsewhere. References:
  //   Dunn P.K., Smyth G.K. (2005). "Series evaluation of Tweedie
  //     exponential dispersion model densities". Statistics and
  //     Computing 15(4):267-280.
  //   Wood S.N. (2017). "Generalized Additive Models" 2nd ed.,
  //     Section 3.5.4.
  //   Jorgensen B. (1997). "The Theory of Dispersion Models",
  //     Chapter 4.

  // Log of the Dunn--Smyth series W(y, phi, p) = sum_{j>=1} W_j,
  // with
  //   W_j = y^(j*alpha_pos) * (p-1)^(-j*alpha_pos)
  //         * phi^(-j*(1+alpha_pos)) * (2-p)^(-j)
  //         * 1 / (j! * Gamma(j*alpha_pos))
  // and alpha_pos := (2-p)/(p-1) > 0 for 1 < p < 2. The summand is
  // unimodal in j; we iterate forward from j = 1, track the running
  // maximum, and break when the term drops more than 37 nats below
  // the maximum (~ 1e-16 relative contribution).
  real tweedie_log_W_series(real y, real phi, real p) {
    real alpha_pos = (2.0 - p) / (p - 1.0);
    real log_z = alpha_pos * log(y)
                 - alpha_pos * log(p - 1.0)
                 - (1.0 + alpha_pos) * log(phi)
                 - log(2.0 - p);
    real log_W = negative_infinity();
    real max_log_W = negative_infinity();
    int max_iter = 1000;
    int passed_peak = 0;
    for (j_loop in 1:max_iter) {
      real j = j_loop;
      real term = j * log_z - lgamma(j + 1.0) - lgamma(j * alpha_pos);
      if (term > max_log_W) {
        max_log_W = term;
      } else {
        passed_peak = 1;
      }
      log_W = log_sum_exp(log_W, term);
      if (passed_peak == 1 && term < max_log_W - 37.0) {
        break;
      }
    }
    return log_W;
  }

  // Exact log-density via the Dunn--Smyth series. Valid for y > 0.
  real tweedie_log_f_series(real y, real mu, real phi, real p) {
    real theta = pow(mu, 1.0 - p) / (1.0 - p);
    real kappa = pow(mu, 2.0 - p) / (2.0 - p);
    return (y * theta - kappa) / phi - log(y)
           + tweedie_log_W_series(y, phi, p);
  }

  // Saddlepoint log-density (Wood 2017, Eq. 3.46 specialized to
  // Tweedie). Valid for y > 0; used outside the series central
  // region for gradient stability when |p - 1.5| >= tau.
  real tweedie_log_f_saddlepoint(real y, real mu, real phi, real p) {
    real deviance = 2.0 * (
      y * (pow(y, 1.0 - p) - pow(mu, 1.0 - p)) / (1.0 - p)
      - (pow(y, 2.0 - p) - pow(mu, 2.0 - p)) / (2.0 - p)
    );
    return -0.5 * log(2.0 * pi() * phi * pow(y, p))
           - deviance / (2.0 * phi);
  }

  // Public log-pdf entry point. The y = 0 case uses the closed-form
  // compound Poisson mass log Pr(Y = 0) = -mu^(2-p) / (phi * (2-p));
  // the y > 0 case dispatches to series or saddlepoint. The function
  // is declared with comma-separated arguments (Stan user-defined
  // _lpdf convention); the canonical sampling syntax
  // `y ~ tweedie(mu, phi, p)` and the explicit
  // `tweedie_lpdf(y | mu, phi, p)` invocation form both dispatch
  // through this signature.
  real tweedie_lpdf(real y, real mu, real phi, real p) {
    real tau = 0.4;
    real lambda = pow(mu, 2.0 - p) / (phi * (2.0 - p));
    if (y == 0.0) {
      return -lambda;
    }
    if (abs(p - 1.5) < tau) {
      return tweedie_log_f_series(y, mu, phi, p);
    } else {
      return tweedie_log_f_saddlepoint(y, mu, phi, p);
    }
  }

  // Compound Poisson-gamma random-number generator.
  // N ~ Poisson(lambda); Y = sum_{i=1}^N G_i with
  // G_i ~ Gamma(shape, rate), shape = (2-p)/(p-1) > 0,
  // rate = 1 / (phi * (p-1) * mu^(p-1)).
  real tweedie_rng(real mu, real phi, real p) {
    real lambda = pow(mu, 2.0 - p) / (phi * (2.0 - p));
    int N = poisson_rng(lambda);
    if (N == 0) {
      return 0.0;
    }
    real shape = (2.0 - p) / (p - 1.0);
    real rate = 1.0 / (phi * (p - 1.0) * pow(mu, p - 1.0));
    return gamma_rng(N * shape, rate);
  }

  // Sub-phase 8.3.7 (2026-05-21): heterogeneous slot inverse-link
  // dispatch. The data field inv_link_id_per_slot[k] (populated R-side
  // by .gdpar_compute_inv_link_id_per_slot under the L1 refined rule)
  // selects the canonical inverse-link of slot k:
  //   0 = identity (e.g. Gaussian mu, Student-t mu, lognormal_loc_scale mu)
  //   1 = inv_logit (e.g. Beta mu, Bernoulli mu)
  //   2 = exp (e.g. Gaussian sigma, NB phi, Beta phi, Gamma mu/shape,
  //            Tweedie mu/phi, lognormal_loc_scale sigma, mixtures mu/phi)
  // Used by the K = 2 likelihood branches (family_id_k[1] in {1, 3, 5,
  // 6, 7}) to compose the auxiliary slot's transformation under the
  // L1 rule. The K = 3+ families (Student-t, Tweedie, mixtures) retain
  // their canonical Stan idioms (*_log_lpmf, bernoulli_logit_lpmf) and
  // a separate R-side guard rejects heterogeneous location families
  // outside {1, 3, 5, 6, 7} until a future sub-phase extends the
  // refactor.
  real apply_inv_link_by_id(int link_id, real eta) {
    if (link_id == 0) return eta;
    if (link_id == 1) return inv_logit(eta);
    return exp(eta);
  }

  // {{CANONICAL_HELPERS}}
}

data {
  int<lower=1> n;
  int<lower=2> K;                       // number of K-individual slots
                                        // (Unit 3 only fires when K > 1)

  array[K] int<lower=1, upper=13> family_id_k;
  // K = 1 templates use a single scalar family_id; here each slot
  // carries its own. Sub-phase 8.3.7 (2026-05-21) admits heterogeneous
  // families per slot for the K = 2 location families {1 = gaussian,
  // 3 = neg_binomial_2, 5 = beta, 6 = gamma, 7 = lognormal_loc_scale};
  // for the K = 3+ families (8 = student_t, 9 = tweedie, 10..13 =
  // mixtures) heterogeneity is deferred to a future sub-phase. The
  // likelihood block still branches on family_id_k[1] (the location
  // slot); the auxiliary slots consume inv_link_id_per_slot[k] via
  // apply_inv_link_by_id() to apply the slot-specific canonical link.

  array[K] int<lower=0, upper=2> inv_link_id_per_slot;
  // R-side helper .gdpar_compute_inv_link_id_per_slot populates this
  // vector under the L1 refined rule (project memory
  // project_gdpar_block_8_3_7_closed, decision D3.5):
  //   - homogeneous slot k (family_id_k[k] == family_id_k[1]): the ID
  //     reflects the canonical link of slot k of the location family;
  //   - heterogeneous slot k (family_id_k[k] != family_id_k[1]): the
  //     ID reflects the canonical link of slot 1 of family_id_k[k].
  // Values: 0 = identity, 1 = inv_logit, 2 = exp. Consumed by the K
  // = 2 likelihood branches (location family in {1, 3, 5, 6, 7}).

  array[K] int<lower=0, upper=1> use_a_k;
  array[K] int<lower=0, upper=1> use_b_k;
  int<lower=0, upper=1> use_W;          // W is globally shared

  array[K] int<lower=0> J_a_per_k;      // J_a^{(k)} per slot
  array[K] int<lower=0> J_b_per_k;      // J_b^{(k)} per slot
  int<lower=0> J_a_max;                 // max over k of J_a_per_k[k]
  int<lower=0> J_b_max;                 // max over k of J_b_per_k[k]

  // Padded per-slot design matrices. Coordinate p = 1 implies a
  // single design matrix per slot; the ragged column counts are
  // handled by the same flat-pack / pad-to-max device used in
  // amm_distrib_multi.stan.
  array[K] matrix[n, J_a_max] Z_a_k;
  array[K] matrix[n, J_b_max] Z_b_k;

  // Global modulating component.
  int<lower=0> dim_W;
  int<lower=0> d;
  matrix[n, d] X;

  // Outcomes. p = 1 implies a single column per slot; the outcome
  // y is the same across slots (the K slots parameterize *one*
  // univariate response). Both real and integer storage are
  // allocated; the consumed branch follows family_id_k[1].
  vector[n] y_real;
  array[n] int y_int;

  // theta_anchor per slot. p = 1 implies one scalar per slot. The
  // R-side packager exposes this as a vector of length K with names
  // matching the slot order of the gdpar_formula_set.
  vector[K] theta_anchor_K;

  // Dispersion flags carried per slot for inter-slot heterogeneity
  // (e.g. Gaussian K = 2 has no sigma_y entry because the scale is
  // K-individual; neg_binomial_2 K = 2 has no phi entry for the
  // same reason). The flags reproduce the K = 1 + p >= 1 behaviour
  // on slots whose param_spec retains population scope, kept here
  // for future extension to mixed (K-individual + population) family
  // patterns.
  array[K] int<lower=0, upper=1> use_dispersion_y_k;
  array[K] int<lower=0, upper=1> use_dispersion_phi_k;
  // Sub-phase 8.3.8 (2026-05-22): see amm_main.stan header.
  int<lower=0, upper=2> W_type_id;
  int<lower=0> W_n_knots_full;
  vector[W_n_knots_full] W_knots_full;
  int<lower=0> W_degree;

  // Grouping (Block 6.5). The grouping structure is shared across
  // slots: a single (J_groups, group_id) tuple drives the
  // hierarchical anchor on theta_ref_k[k][g] for every k.
  int<lower=0, upper=1> use_groups;
  int<lower=1> J_groups;
  array[n] int<lower=1, upper=J_groups> group_id;{{DATA_CP_A_PER_K_DECL}}

  // Sub-phase 8.6.B / 8.6.C convention: every EB template declares
  // K_slots and p_dim explicitly so the R-side .gdpar_eb_* helpers can
  // populate them uniformly. Under Path B K_slots == K and p_dim == 1
  // (univariate population reference per slot); the body of this
  // template does not branch on either field.
  int<lower=1> K_slots;
  int<lower=1> p_dim;
}

transformed data {
  array[K] int J_a_free;
  array[K] int J_b_free;
  array[K] int a_raw_offset;
  array[K] int c_b_raw_offset;
  int total_J_a_free = 0;
  int total_J_b_free = 0;
  int dim_W_eff = (use_W == 1) ? dim_W : 0;
  int d_eff = (use_W == 1) ? d : 0;
  // sigma_a / sigma_b carry one entry per active slot. Following
  // the K = 1 + p >= 1 template, the sizes collapse to zero when
  // the component is globally off (use_a_k all zeros — captured
  // here by the conjunctive flag any_use_a).
  int any_use_a = 0;
  int any_use_b = 0;
  int sigma_W_size;
  // family-driven dispersion sizing per slot. Each slot may carry
  // sigma_y (Gaussian, location only) or phi (NB, location only)
  // depending on whether the scale parameter of the slot is
  // population-scoped (handled here) or K-individual (modeled via
  // its own AMM in another slot).
  int sigma_y_size = 0;
  int phi_size = 0;

  for (k in 1:K) {
    J_a_free[k] = (use_a_k[k] == 1 && J_a_per_k[k] > 0) ? J_a_per_k[k] : 0;
    J_b_free[k] = (use_b_k[k] == 1 && J_b_per_k[k] > 0) ? J_b_per_k[k] : 0;
    a_raw_offset[k] = total_J_a_free;
    c_b_raw_offset[k] = total_J_b_free;
    total_J_a_free += J_a_free[k];
    total_J_b_free += J_b_free[k];
    if (use_a_k[k] == 1) any_use_a = 1;
    if (use_b_k[k] == 1) any_use_b = 1;
    sigma_y_size += use_dispersion_y_k[k];
    phi_size += use_dispersion_phi_k[k];
  }
  sigma_W_size = (use_W == 1 && dim_W > 0) ? 1 : 0;
  // Option A (RG.6, D96): compact sigma_a_k to slots carrying free a
  // coefficients (J_a_free > 0). A slot with no free a coefficient ---
  // either no a() at all (use_a_k == 0, e.g. phi ~ 1) or an
  // intercept-only a() (a = ~ 1) fully absorbed into theta_ref ---
  // would otherwise declare a sampled-but-unused half-prior scale: a
  // flat direction (non-identification) that stalls the sampler.
  // sigma_a_idx maps each slot to its position in the compacted vector
  // (0 if the slot carries no free a coefficient). When every slot
  // carries free a coefficients, n_sigma_a == K and sigma_a_idx is the
  // identity, so the generated model is mathematically identical to the
  // pre-compaction one (the goldens stay bit-identical).
  int n_sigma_a = 0;
  array[K] int sigma_a_idx;
  for (k in 1:K) {
    if (J_a_free[k] > 0) {
      n_sigma_a += 1;
      sigma_a_idx[k] = n_sigma_a;
    } else {
      sigma_a_idx[k] = 0;
    }
  }
}

parameters {
  // Per-slot anchors. p = 1 implies a scalar anchor per slot per
  // group. The container is array[J_groups] vector[K] to mirror
  // the array[J_groups] vector[p] declaration of the K = 1 + p >= 1
  // template, which keeps the model block code structurally close
  // across the two regimes.
  array[J_groups] vector[K] theta_ref_k;

  // Hierarchical hyperparameters of the per-group anchor prior.
  // Coord-wise (here slot-wise) i.i.d. when use_groups == 1.
  array[use_groups] vector[K] mu_theta_ref_k;
  array[use_groups] vector<lower=0>[K] sigma_theta_ref_k;

  // Per-slot scale of the additive component (compacted to slots with
  // free a coefficients; see sigma_a_idx in transformed data, D96).
  array[n_sigma_a] real<lower=0> sigma_a_k;
  vector[total_J_a_free] a_raw;

  // Per-slot scale of the multiplicative component (linear
  // reparametrization on c_b_k = theta_ref_k * b_coef_k).
  array[K * any_use_b] real<lower=0> sigma_b_k;
  vector[total_J_b_free] c_b_k_raw;

  // Global modulating component.
  array[sigma_W_size] real<lower=0> sigma_W;
  matrix[dim_W_eff, d_eff] W_raw;

  // Population-scoped dispersion entries (only emitted for slots
  // whose param_spec.scope remains population after promotion).
  array[sigma_y_size] real<lower=0> sigma_y_pop;
  array[phi_size] real<lower=0> phi_pop;
}

transformed parameters {
  // Reconstruct per-slot a_coef_k / c_b_k / b_coef_k from the
  // flat-packed raw vectors using the per-slot offsets and the
  // active flag any_use_a / any_use_b.
  array[K] vector[J_a_max] a_coef_k;
  array[K] vector[J_b_max] c_b_k;
  array[K] vector[J_b_max] b_coef_k;
  matrix[n, K] eta_k;

  for (k in 1:K) {
    a_coef_k[k] = rep_vector(0, J_a_max);
    c_b_k[k] = rep_vector(0, J_b_max);
    b_coef_k[k] = rep_vector(0, J_b_max);
  }

{{TP_A_BLOCK}}

  if (any_use_b == 1) {
    for (k in 1:K) {
      if (use_b_k[k] == 1 && J_b_per_k[k] > 0) {
        for (j in 1:J_b_free[k]) {
          c_b_k[k][j] = c_b_k_raw[c_b_raw_offset[k] + j] * sigma_b_k[k];
        }
        // b_coef_k[k] = c_b_k[k] / theta_ref_k[1][k] is reported as a
        // derived quantity only in the single-anchor regime. With
        // per-group anchors b_coef_k stays at zero; downstream callers
        // derive it per group from c_b_k and theta_ref_k[g][k] as in
        // the K = 1 + p >= 1 template.
        if (use_groups == 0 && theta_ref_k[1][k] != 0) {
          for (j in 1:J_b_per_k[k]) {
            b_coef_k[k][j] = c_b_k[k][j] / theta_ref_k[1][k];
          }
        }
      }
    }
  }

  for (i in 1:n) {
    int g_i = group_id[i];
    for (k in 1:K) {
      real theta_ref_ik = theta_ref_k[g_i][k];
      real eta_ik = theta_ref_ik;
      if (use_a_k[k] == 1 && J_a_per_k[k] > 0) {
        eta_ik += dot_product(
          to_vector(Z_a_k[k][i, 1:J_a_per_k[k]]),
          a_coef_k[k][1:J_a_per_k[k]]
        );
      }
      if (use_b_k[k] == 1 && J_b_per_k[k] > 0) {
        eta_ik += dot_product(
          to_vector(Z_b_k[k][i, 1:J_b_per_k[k]]),
          c_b_k[k][1:J_b_per_k[k]]
        );
      }
      if (use_W == 1 && dim_W > 0 && d > 0) {
        // W is globally shared across slots. Each slot k draws
        // its modulating contribution from the *same* basis
        // evaluated at its own theta_ref_k and anchored at its
        // own theta_anchor_K[k]. The basis dimension dim_W refers
        // to the number of basis functions; for K > 1 + p = 1 the
        // basis is univariate (per-slot scalar input), so the row
        // count of W_raw is dim_W rather than the multivariate
        // p-stride dim_W * p used in amm_distrib_multi.stan.
        // Sub-phase 8.3.8 (2026-05-22): basis dispatch via
        // apply_W_basis_diff() supports W_type_id in {1=polynomial,
        // 2=bspline} uniformly across the K slots.
        vector[d] W_diff_x_k = rep_vector(0, d);
        vector[dim_W] basis_diff_k = apply_W_basis_diff(
          W_type_id, theta_ref_ik, theta_anchor_K[k], dim_W,
          W_degree, W_n_knots_full, W_knots_full
        );
        for (jj_local in 1:dim_W) {
          for (jjj in 1:d) {
            W_diff_x_k[jjj] += basis_diff_k[jj_local] * W_raw[jj_local, jjj]{{W_SCALE}};
          }
        }
        eta_ik += dot_product(W_diff_x_k, to_vector(X[i]));
      }
      eta_k[i, k] = eta_ik;
    }
  }
}

model {
  // Anchor priors per slot. Hierarchical when use_groups == 1;
  // flat single-anchor i.i.d. across slots otherwise. The block is
  // emitted by R/stan_codegen.R via the THETA_REF_PRIOR_BLOCK
  // placeholder so families with structurally tight slots (e.g.,
  // Tweedie's p in (1.01, 1.99), sub-phase 8.3.5b) can apply a
  // slot-specific prior (E7 of D5: slice notation + uniform on the
  // power slot) without giving up the vectorized canonical prior on
  // the remaining slots.
  {{THETA_REF_PRIOR_BLOCK}}

  // Per-slot scale priors.
  if (any_use_a == 1) {
    sigma_a_k ~ {{PRIOR_SIGMA_A}};
{{MODEL_A_BLOCK}}
  }
  if (any_use_b == 1) {
    sigma_b_k ~ {{PRIOR_SIGMA_B}};
    c_b_k_raw ~ normal(0, 1);
  }
  if (use_W == 1 && dim_W > 0) {
    sigma_W[1] ~ {{PRIOR_SIGMA_W}};
    to_vector(W_raw) ~ {{W_PRIOR}};
  }
  if (sigma_y_size > 0) {
    sigma_y_pop ~ {{PRIOR_SIGMA_Y}};
  }
  if (phi_size > 0) {
    phi_pop ~ {{PRIOR_PHI}};
  }

  // Distributional regression likelihood. The K = 2 case is
  // wired explicitly per built-in family. Slot 1 is the location
  // slot (link applied via family-specific inverse); slot 2 is the
  // dispersion slot (link applied via exp()). The homogeneous-family
  // assumption (family_id_k all equal) holds in Unit 3; family_id_k[1]
  // drives the dispatch.
  if (family_id_k[1] == 1) {
    // Gaussian K = 2: y ~ Normal(mu, sigma) where mu = eta_k[i, 1]
    // (slot 1 identity, canonical for Gaussian location) and sigma is
    // produced by apply_inv_link_by_id(inv_link_id_per_slot[2],
    // eta_k[i, 2]). Under homogeneity (family_id_k[2] == 1) the slot
    // 2 link is log per Gaussian's canonical sigma parametrization
    // (inv_link_id_per_slot[2] == 2 -> exp). Under heterogeneity
    // (family_id_k[2] in {5 beta}) the slot 2 link follows L1 (slot 1
    // of the heterogeneous family); e.g. Beta yields inv_logit and
    // sigma in (0, 1).
    for (i in 1:n) {
      real sigma_i = apply_inv_link_by_id(inv_link_id_per_slot[2],
                                          eta_k[i, 2]);
      y_real[i] ~ normal(eta_k[i, 1], sigma_i);
    }
  } else if (family_id_k[1] == 3) {
    // neg_binomial_2 K = 2: y_int ~ NegBinomial2(mu, phi) with mu =
    // exp(eta_k[i, 1]) (Stan idiom via neg_binomial_2 form) and phi =
    // apply_inv_link_by_id(inv_link_id_per_slot[2], eta_k[i, 2]). Under
    // homogeneity slot 2 = phi log (inv_link_id_per_slot[2] == 2 ->
    // exp), reproducing the legacy neg_binomial_2_log(eta_mu,
    // exp(eta_phi)) up to factoring of the eta_mu exponentiation. Under
    // heterogeneity (e.g. Beta in slot 2) phi follows the
    // heterogeneous family's slot 1 link.
    for (i in 1:n) {
      real phi_i = apply_inv_link_by_id(inv_link_id_per_slot[2],
                                        eta_k[i, 2]);
      y_int[i] ~ neg_binomial_2(exp(eta_k[i, 1]), phi_i);
    }
  } else if (family_id_k[1] == 5) {
    // Beta K = 2: y ~ BetaProportion(inv_logit(eta_mu), phi).
    // Canonical mean-precision parametrization (mu logit, phi log).
    // Under homogeneity (family_id_k[2] == 5) phi = exp(eta_k[i, 2]).
    // Under heterogeneity phi follows the heterogeneous family's slot
    // 1 inverse-link.
    for (i in 1:n) {
      real phi_i = apply_inv_link_by_id(inv_link_id_per_slot[2],
                                        eta_k[i, 2]);
      y_real[i] ~ beta_proportion(inv_logit(eta_k[i, 1]), phi_i);
    }
  } else if (family_id_k[1] == 6) {
    // Gamma K = 2: y ~ Gamma(shape, shape / mu). Canonical
    // mean-shape parametrization (mu log, shape log); mu =
    // exp(eta_k[i, 1]) is fixed (Gamma location), shape =
    // apply_inv_link_by_id(inv_link_id_per_slot[2], eta_k[i, 2]).
    // Stan's gamma uses (alpha, beta) where mean = alpha / beta; we
    // convert via beta = alpha / mu.
    for (i in 1:n) {
      real mu_i    = exp(eta_k[i, 1]);
      real shape_i = apply_inv_link_by_id(inv_link_id_per_slot[2],
                                          eta_k[i, 2]);
      y_real[i] ~ gamma(shape_i, shape_i / mu_i);
    }
  } else if (family_id_k[1] == 7) {
    // Lognormal_loc_scale K = 2 via descriptor-based custom-family
    // wiring (D-A3.B, sub-phase 8.3.4). Canonical (mu identity, sigma
    // log) on the log-scale of y. sigma =
    // apply_inv_link_by_id(inv_link_id_per_slot[2], eta_k[i, 2]).
    for (i in 1:n) {
      real sigma_i = apply_inv_link_by_id(inv_link_id_per_slot[2],
                                          eta_k[i, 2]);
      y_real[i] ~ lognormal(eta_k[i, 1], sigma_i);
    }
  } else if (family_id_k[1] == 8) {
    // Student-t K = 3 (sub-phase 8.3.5a). Canonical (mu, log sigma,
    // log nu) parametrization: slot 1 carries the location on the
    // identity link, slot 2 carries log(sigma), slot 3 carries
    // log(nu). Stan's student_t(nu, mu, sigma) consumes the raw
    // degrees-of-freedom (positive real) as its first argument.
    for (i in 1:n) {
      y_real[i] ~ student_t(exp(eta_k[i, 3]),
                            eta_k[i, 1],
                            exp(eta_k[i, 2]));
    }
  } else if (family_id_k[1] == 9) {
    // Tweedie K = 3 (sub-phase 8.3.5b). Canonical (mu, log phi, p)
    // parametrization: slot 1 carries log(mu) (mean on the log link),
    // slot 2 carries log(phi) (dispersion on the log link), slot 3
    // carries p (power on the identity link, structurally bounded to
    // (1.01, 1.99) by the slot-specific prior emitted by R-side per
    // decision E7 of D5). The custom tweedie_lpdf defined above
    // dispatches Dunn--Smyth series or saddlepoint depending on the
    // hybrid boundary |p - 1.5| < 0.4.
    for (i in 1:n) {
      y_real[i] ~ tweedie(exp(eta_k[i, 1]),
                          exp(eta_k[i, 2]),
                          eta_k[i, 3]);
    }
  } else if (family_id_k[1] == 10) {
    // ZIP K = 2 (sub-phase 8.3.6). Canonical (mu log, pi logit)
    // parametrization: slot 1 carries log(mu) of the Poisson count
    // component, slot 2 carries logit(pi) of the structural zero
    // probability. Likelihood (Lambert 1992):
    //   P(y = 0) = pi + (1 - pi) * Poisson(0 | mu)
    //   P(y > 0) = (1 - pi) * Poisson(y | mu)
    // Stable Stan implementation: log_sum_exp on the y = 0 branch
    // (structural zero vs Poisson sampling zero) and
    // bernoulli_logit_lpmf(0 | eta_pi) + poisson_log_lpmf(y | eta_mu)
    // on the y > 0 branch.
    for (i in 1:n) {
      real eta_mu_i = eta_k[i, 1];
      real eta_pi_i = eta_k[i, 2];
      if (y_int[i] == 0) {
        target += log_sum_exp(
          bernoulli_logit_lpmf(1 | eta_pi_i),
          bernoulli_logit_lpmf(0 | eta_pi_i)
            + poisson_log_lpmf(0 | eta_mu_i)
        );
      } else {
        target += bernoulli_logit_lpmf(0 | eta_pi_i)
                + poisson_log_lpmf(y_int[i] | eta_mu_i);
      }
    }
  } else if (family_id_k[1] == 11) {
    // ZINB K = 3 (sub-phase 8.3.6). Canonical (mu log, phi log, pi
    // logit) parametrization: slot 1 = log(mu) of the NB count, slot
    // 2 = log(phi) of the NB dispersion, slot 3 = logit(pi) of the
    // structural zero probability. Likelihood (Greene 1994):
    //   P(y = 0) = pi + (1 - pi) * NegBinomial2(0 | mu, phi)
    //   P(y > 0) = (1 - pi) * NegBinomial2(y | mu, phi)
    for (i in 1:n) {
      real eta_mu_i  = eta_k[i, 1];
      real phi_i     = exp(eta_k[i, 2]);
      real eta_pi_i  = eta_k[i, 3];
      if (y_int[i] == 0) {
        target += log_sum_exp(
          bernoulli_logit_lpmf(1 | eta_pi_i),
          bernoulli_logit_lpmf(0 | eta_pi_i)
            + neg_binomial_2_log_lpmf(0 | eta_mu_i, phi_i)
        );
      } else {
        target += bernoulli_logit_lpmf(0 | eta_pi_i)
                + neg_binomial_2_log_lpmf(y_int[i] | eta_mu_i, phi_i);
      }
    }
  } else if (family_id_k[1] == 12) {
    // Hurdle-Poisson K = 2 (sub-phase 8.3.6). Canonical (mu log, pi
    // logit) parametrization. Distinct from ZIP in that the zero is
    // a structural decision (Bernoulli hurdle) and the positive
    // branch is a Poisson truncated at one (Mullahy 1986):
    //   P(y = 0)  = pi
    //   P(y > 0)  = (1 - pi) * Poisson(y | mu) / (1 - exp(-mu))
    // log1m_exp(-exp(eta_mu)) computes log(1 - Poisson(0 | mu)) in a
    // numerically stable way for the truncation denominator.
    for (i in 1:n) {
      real eta_mu_i = eta_k[i, 1];
      real eta_pi_i = eta_k[i, 2];
      if (y_int[i] == 0) {
        target += bernoulli_logit_lpmf(1 | eta_pi_i);
      } else {
        target += bernoulli_logit_lpmf(0 | eta_pi_i)
                + poisson_log_lpmf(y_int[i] | eta_mu_i)
                - log1m_exp(-exp(eta_mu_i));
      }
    }
  } else if (family_id_k[1] == 13) {
    // Hurdle-NB K = 3 (sub-phase 8.3.6). Canonical (mu log, phi log,
    // pi logit) parametrization. Symmetric to Hurdle-Poisson with
    // NegBinomial2 in the positive branch truncated at one:
    //   P(y = 0) = pi
    //   P(y > 0) = (1 - pi) * NegBinomial2(y | mu, phi) /
    //              (1 - NegBinomial2(0 | mu, phi))
    for (i in 1:n) {
      real eta_mu_i = eta_k[i, 1];
      real phi_i    = exp(eta_k[i, 2]);
      real eta_pi_i = eta_k[i, 3];
      if (y_int[i] == 0) {
        target += bernoulli_logit_lpmf(1 | eta_pi_i);
      } else {
        target += bernoulli_logit_lpmf(0 | eta_pi_i)
                + neg_binomial_2_log_lpmf(y_int[i] | eta_mu_i, phi_i)
                - log1m_exp(neg_binomial_2_log_lpmf(0 | eta_mu_i,
                                                    phi_i));
      }
    }
  }
  // family_id_k[1] in {2 (poisson), 4 (bernoulli)}: these are
  // K = 1 families (single eligible param_spec); the dispatcher
  // never routes them here. The guard in R/families.R ensures
  // promotion to per_observation requires the slot name to be
  // eligible.
}

generated quantities {
  // Per-observation log-likelihood, fitted linear predictors per
  // slot, and posterior-predictive draws. The shape mirrors the
  // K = 1 + p >= 1 template up to the swap of p with K and the
  // distributional-regression likelihood selection.
  vector[n] log_lik;
  matrix[n, K] theta_i_k = eta_k;
  vector[n] y_pred;

  for (i in 1:n) {
    if (family_id_k[1] == 1) {
      real sigma_i = apply_inv_link_by_id(inv_link_id_per_slot[2],
                                          eta_k[i, 2]);
      log_lik[i] = normal_lpdf(y_real[i] | eta_k[i, 1], sigma_i);
      y_pred[i] = normal_rng(eta_k[i, 1], sigma_i);
    } else if (family_id_k[1] == 3) {
      real mu_i  = exp(eta_k[i, 1]);
      real phi_i = apply_inv_link_by_id(inv_link_id_per_slot[2],
                                        eta_k[i, 2]);
      log_lik[i] = neg_binomial_2_lpmf(y_int[i] | mu_i, phi_i);
      y_pred[i]  = neg_binomial_2_rng(mu_i, phi_i);
    } else if (family_id_k[1] == 5) {
      real mu_i  = inv_logit(eta_k[i, 1]);
      real phi_i = apply_inv_link_by_id(inv_link_id_per_slot[2],
                                        eta_k[i, 2]);
      log_lik[i] = beta_proportion_lpdf(y_real[i] | mu_i, phi_i);
      y_pred[i]  = beta_proportion_rng(mu_i, phi_i);
    } else if (family_id_k[1] == 6) {
      real mu_i    = exp(eta_k[i, 1]);
      real shape_i = apply_inv_link_by_id(inv_link_id_per_slot[2],
                                          eta_k[i, 2]);
      log_lik[i] = gamma_lpdf(y_real[i] | shape_i, shape_i / mu_i);
      y_pred[i]  = gamma_rng(shape_i, shape_i / mu_i);
    } else if (family_id_k[1] == 7) {
      real sigma_i = apply_inv_link_by_id(inv_link_id_per_slot[2],
                                          eta_k[i, 2]);
      log_lik[i] = lognormal_lpdf(y_real[i] | eta_k[i, 1], sigma_i);
      y_pred[i]  = lognormal_rng(eta_k[i, 1], sigma_i);
    } else if (family_id_k[1] == 8) {
      real nu_i    = exp(eta_k[i, 3]);
      real mu_i    = eta_k[i, 1];
      real sigma_i = exp(eta_k[i, 2]);
      log_lik[i] = student_t_lpdf(y_real[i] | nu_i, mu_i, sigma_i);
      y_pred[i] = student_t_rng(nu_i, mu_i, sigma_i);
    } else if (family_id_k[1] == 9) {
      real mu_i  = exp(eta_k[i, 1]);
      real phi_i = exp(eta_k[i, 2]);
      real p_i   = eta_k[i, 3];
      log_lik[i] = tweedie_lpdf(y_real[i] | mu_i, phi_i, p_i);
      y_pred[i] = tweedie_rng(mu_i, phi_i, p_i);
    } else if (family_id_k[1] == 10) {
      // ZIP K = 2: log_lik mirrors the model-block dispatch; y_pred
      // draws the structural zero indicator first and then a Poisson
      // sample for the non-zero branch (Lambert 1992).
      real eta_mu_i = eta_k[i, 1];
      real eta_pi_i = eta_k[i, 2];
      if (y_int[i] == 0) {
        log_lik[i] = log_sum_exp(
          bernoulli_logit_lpmf(1 | eta_pi_i),
          bernoulli_logit_lpmf(0 | eta_pi_i)
            + poisson_log_lpmf(0 | eta_mu_i)
        );
      } else {
        log_lik[i] = bernoulli_logit_lpmf(0 | eta_pi_i)
                   + poisson_log_lpmf(y_int[i] | eta_mu_i);
      }
      if (bernoulli_logit_rng(eta_pi_i) == 1) {
        y_pred[i] = 0;
      } else {
        y_pred[i] = poisson_log_rng(eta_mu_i);
      }
    } else if (family_id_k[1] == 11) {
      // ZINB K = 3: symmetric to ZIP with neg_binomial_2_log.
      real eta_mu_i = eta_k[i, 1];
      real phi_i    = exp(eta_k[i, 2]);
      real eta_pi_i = eta_k[i, 3];
      if (y_int[i] == 0) {
        log_lik[i] = log_sum_exp(
          bernoulli_logit_lpmf(1 | eta_pi_i),
          bernoulli_logit_lpmf(0 | eta_pi_i)
            + neg_binomial_2_log_lpmf(0 | eta_mu_i, phi_i)
        );
      } else {
        log_lik[i] = bernoulli_logit_lpmf(0 | eta_pi_i)
                   + neg_binomial_2_log_lpmf(y_int[i] | eta_mu_i,
                                              phi_i);
      }
      if (bernoulli_logit_rng(eta_pi_i) == 1) {
        y_pred[i] = 0;
      } else {
        y_pred[i] = neg_binomial_2_log_rng(eta_mu_i, phi_i);
      }
    } else if (family_id_k[1] == 12) {
      // Hurdle-Poisson K = 2: log_lik mirrors the model-block
      // dispatch. y_pred uses rejection sampling on Poisson(mu) to
      // realize the truncated-at-one positive branch (Mullahy 1986);
      // bounded at 10000 iterations for numerical safety when mu is
      // pathologically small (in which case the family is borderline
      // identifiable in the first place).
      real eta_mu_i = eta_k[i, 1];
      real eta_pi_i = eta_k[i, 2];
      if (y_int[i] == 0) {
        log_lik[i] = bernoulli_logit_lpmf(1 | eta_pi_i);
      } else {
        log_lik[i] = bernoulli_logit_lpmf(0 | eta_pi_i)
                   + poisson_log_lpmf(y_int[i] | eta_mu_i)
                   - log1m_exp(-exp(eta_mu_i));
      }
      if (bernoulli_logit_rng(eta_pi_i) == 1) {
        y_pred[i] = 0;
      } else {
        int y_draw = 0;
        int iter = 0;
        while (y_draw == 0 && iter < 10000) {
          y_draw = poisson_log_rng(eta_mu_i);
          iter = iter + 1;
        }
        y_pred[i] = y_draw;
      }
    } else if (family_id_k[1] == 13) {
      // Hurdle-NB K = 3: symmetric to Hurdle-Poisson with
      // neg_binomial_2_log on the truncated positive branch.
      real eta_mu_i = eta_k[i, 1];
      real phi_i    = exp(eta_k[i, 2]);
      real eta_pi_i = eta_k[i, 3];
      if (y_int[i] == 0) {
        log_lik[i] = bernoulli_logit_lpmf(1 | eta_pi_i);
      } else {
        log_lik[i] = bernoulli_logit_lpmf(0 | eta_pi_i)
                   + neg_binomial_2_log_lpmf(y_int[i] | eta_mu_i,
                                              phi_i)
                   - log1m_exp(neg_binomial_2_log_lpmf(0 | eta_mu_i,
                                                       phi_i));
      }
      if (bernoulli_logit_rng(eta_pi_i) == 1) {
        y_pred[i] = 0;
      } else {
        int y_draw = 0;
        int iter = 0;
        while (y_draw == 0 && iter < 10000) {
          y_draw = neg_binomial_2_log_rng(eta_mu_i, phi_i);
          iter = iter + 1;
        }
        y_pred[i] = y_draw;
      }
    } else {
      log_lik[i] = negative_infinity();
      y_pred[i] = 0;
    }
  }
}
