R中向量化S4类：卷积似然计算的提速方法求助

我正在计算对数似然，该似然是两个分布的卷积，且分布参数取决于对应数据点的值。更准确地说，数据服从的似然是两个Gamma分布的和（注意：这里是随机变量的算术运算，不是混合模型），二者具有独立的形状参数和尺度参数，且形状参数是数据的函数，数学表达式为：
$y_i = p\text{Gamma}(k_1+b_1x_i, t_1) + (1-p)\text{Gamma}(k_2+b_2x_i, t_2)$

目前使用distr包计算似然的代码速度极慢：

sapply(1:10000, function(i){distr::d(p*Gammad(shape = k1+b1*x[i], scale = t1)+(1-p)*Gammad(shape = k2+b2*x[i], scale = t2))(y[i])})

尝试向量化处理时出现错误：

distr::d(Vectorize(p*Gammad(shape = k1+b1*x, scale = t1)+(1-p)*Gammad(shape = k2+b2*x, scale = t2)))(y)

错误信息：Error in .multm(e1, e2, "AbscontDistribution") : length of operator must be 1

想知道是否可以实现向量化，或者有其他加速计算的方法。

• 放弃distr包的向量化尝试：distr的分布对象不支持批量运算，Vectorize无法处理分布对象的算术操作，因此会触发报错。
• 用数值卷积替代符号分布运算：直接对每个数据点对应的两个Gamma分布概率密度做数值卷积，再插值得到目标点的密度值，避免distr包的低效符号运算。示例代码：

# 定义单个数据点的似然计算函数
single_likelihood <- function(x_i, y_i, p, k1, b1, t1, k2, b2, t2) {
  # 计算当前数据点对应的Gamma形状参数
  shape1 <- k1 + b1 * x_i
  shape2 <- k2 + b2 * x_i
  
  # 确定采样范围（基于y_i和Gamma分布的均值，保证覆盖有效区间）
  mean1 <- shape1 * t1 * p
  mean2 <- shape2 * t2 * (1 - p)
  max_val <- y_i + 3 * sqrt(mean1^2/shape1 + mean2^2/shape2)  # 3倍标准差覆盖
  x_grid <- seq(0.01, max_val, length.out = 200)
  
  # 生成缩放后的Gamma密度（p*Gamma和(1-p)*Gamma的密度变换）
  d1 <- dgamma(x_grid / p, shape = shape1, scale = t1) / p
  d2 <- dgamma(x_grid / (1 - p), shape = shape2, scale = t2) / (1 - p)
  
  # 计算两个密度的卷积
  conv_result <- convolve(d1, rev(d2), type = "open")
  
  # 插值得到y_i对应的密度值
  conv_grid <- seq(x_grid[1] + x_grid[1], x_grid[length(x_grid)] + x_grid[length(x_grid)], length.out = length(conv_result))
  approx(conv_grid, conv_result, xout = y_i)$y
}

# 向量化批量计算似然
likelihoods <- mapply(single_likelihood, x_i = x, y_i = y, 
                      MoreArgs = list(p = p, k1 = k1, b1 = b1, t1 = t1, k2 = k2, b2 = b2, t2 = t2))

• 并行运算进一步加速：针对大样本量，使用多核心并行计算替代单线程的mapply：

library(parallel)
# 创建并行集群
cl <- makeCluster(detectCores() - 1)
# 导出所需变量和函数到集群节点
clusterExport(cl, c("p", "k1", "b1", "t1", "k2", "b2", "t2", "single_likelihood"))
# 并行计算似然
likelihoods <- parSapply(cl, 1:length(x), function(i) {
  single_likelihood(x[i], y[i], p, k1, b1, t1, k2, b2, t2)
})
# 关闭集群
stopCluster(cl)

• 精度与速度的平衡调整：可以根据实际需求调整x_grid的长度（比如减少到150）和采样范围，在保证似然计算精度的前提下降低计算量。

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