如何在rethinking包的ulam接口中实现预测矩阵与参数向量相乘

Got it, I've run into this exact issue when switching from quap to ulam for spline models in Statistical Rethinking too! The key is adjusting how you handle the matrix-vector multiplication since ulam has slightly different syntax requirements for linear algebra operations compared to quap.

Here's the working solution using ulam: #

First, we'll reuse your data prep code, then modify the model for HMC:

data("cherry_blossoms")
d = cherry_blossoms
d2 = d[ complete.cases(d$temp), ]
num_knots = 15
knot_list = quantile(d2$year, probs = seq(0,1,length.out = num_knots))
library(splines)
B = bs(d2$year,knots = knot_list[-c(1,num_knots)],degree = 3, intercept = TRUE) # 1124 x 17 matrix

library(rethinking)

# Option 1: Use rows_dot_product (most intuitive for per-observation calculations)
m4.7_ulam = ulam(
  alist(
    T ~ dnorm(mu, sigma),
    mu <- a + rows_dot_product(B, w),  # Replaces B %*% w with row-wise dot product
    a ~ dnorm(6, 10),
    w ~ dnorm(0, 1),
    sigma ~ dexp(1)
  ), 
  data = list(T = d2$temp, B = B),
  start = list(w = rep(0, ncol(B))),
  chains = 4,
  cores = 4,  # Use multiple cores for faster sampling
  iter = 2000
)

# Option 2: Explicitly declare vector dimension for direct matrix multiplication
K <- ncol(B)
m4.7_ulam_v2 = ulam(
  alist(
    T ~ dnorm(mu, sigma),
    mu <- a + B %*% w,  # Same syntax as quap, but we need to define w's length
    a ~ dnorm(6, 10),
    w[K] ~ dnorm(0, 1),  # Tell ulam w is a vector of length K
    sigma ~ dexp(1)
  ), 
  data = list(T = d2$temp, B = B, K = K),
  start = list(w = rep(0, K)),
  chains = 4,
  cores = 4,
  iter = 2000
)

Why this works: #

• Option 1: rows_dot_product(B, w) is a helper function in the rethinking package designed specifically for this scenario. It calculates the dot product of each row in matrix B with the parameter vector w, which is exactly what B %*% w does—but it's formatted to play nicely with ulam's internal parameter handling.
• Option 2: If you prefer sticking close to your original quap syntax, explicitly declaring the length of w (using w[K]) tells ulam that w is a vector of fixed length, allowing it to correctly interpret the matrix multiplication B %*% w.

Post-Sampling Checks: #

After running the model, make sure to verify convergence:

# Check traceplots to ensure chains mix well
traceplot(m4.7_ulam)

# Summarize posterior parameters (compare to quap results to confirm consistency)
precis(m4.7_ulam, depth = 2)

You should see results that align closely with your original quap model—HMC will converge to the same posterior mode here, just with full posterior samples instead of a quadratic approximation.

内容的提问来源于stack exchange，提问作者Dedula33