Module bayesian 

Bayesian function-on-scalar regression via conjugate Gibbs sampler (REG-06-04).

Fits Y_i(t) = μ(t) + Σ_j x̃_ij · β_j(t) + ε_i(t) where x̃ are the mean-centered scalar predictors and β_j(t) are functional coefficients. The functional dimension is compressed with FPCA (fdata_to_pc_1d): the response is projected onto its top-K functional principal components, and for each component k the FPC scores are regressed on the predictors with a conjugate Normal / Inverse-Gamma Gibbs sampler. Each retained draw reconstructs the coefficient functions β_j(t) = Σ_k b_{jk} · φ_k(t) from the FPCA rotation, so posterior summaries (mean + pointwise credible bands) are obtained directly in the functional domain.

Full conditionals (per FPC component k)

Model: ξ_k = X̃ · b_k + ε_k, ξ_k ∈ Rⁿ (scores of component k), X̃ ∈ R^{n×p}. Priors: b_k | σ²_k ~ N(0, τ²·I_p), σ²_k ~ IG(a₀, b₀).

• b_k | · ~ N(μ_post, A⁻¹) with precision A = X̃'X̃/σ²_k + I_p/τ², μ_post = A⁻¹ · X̃'ξ_k / σ²_k. Sampled via Cholesky A = LLᵀ: b_k = μ_post + Lᵀ⁻¹ z, z ~ N(0, I_p) (Rue 2001) — Cov = (LLᵀ)⁻¹ = A⁻¹.
• σ²_k | · ~ IG(a₀ + n/2, b₀ + RSS_k/2), drawn as 1 / Gamma(α, 1/β).

The chain is fully deterministic given config.seed (StdRng::seed_from_u64(seed)).

References

Rue (2001). Fast sampling of Gaussian Markov random fields. JRSS-B 63(2). Goldsmith et al. (2015). JCGS 23(1). Jiang et al. (2025), arXiv:2505.05633.

Divergences from refund

refund’s Bayesian FOSR uses spline basis priors with random effects; this implementation uses FPCA score compression (fdata_to_pc_1d) for simplicity and zero new dependencies. Pointwise credible bands only (no simultaneous bands).

Functions

bayesian_fosr

Bayesian function-on-scalar regression via conjugate Gibbs sampler.