如何从联合高斯分布定义计算边缘密度函数？R代码报错排查

1. integrate报错原因及修复 #

报错的核心原因是：integrate函数执行时会向被积函数传入向量形式的x（而非单独标量），但你编写的myfun用c(x, y)组合输入，当x是向量时，会生成长度为length(x)+1的一维向量，而dmvnorm要求输入的样本点维度（此处为2维）与均值mu的维度（长度2）一致，因此触发维度不匹配错误。

修复方法是修改myfun，将x和y组合成每行一个二维样本点的矩阵，dmvnorm支持矩阵输入，可自动处理向量形式的x：

library(mvtnorm)

myfun = function(x, y, mu, ht) {
  # 把向量x和重复后的y绑定成n行2列的矩阵
  points <- cbind(x, rep(y, length(x)))
  dmvnorm(points, mean = mu, sigma = ht)
}

现在重新运行你的integrate代码即可正常执行：

integrate(myfun,
          lower = -10,
          upper = 2,
          y = -2,
          mu = c(0.06, 0.03),
          ht = diag(2))

2. 联合高斯边缘密度的计算方法 #

联合高斯分布的边缘密度有两种计算方式：

数值积分法（基于定义）

边缘密度的定义是联合密度对另一变量的积分，比如计算X的边缘密度，需对Y的所有可能值积分联合密度f(x,y)。示例代码如下：

# 定义被积函数：固定x，对y积分
myfun_edge_x = function(y, x, mu, ht) {
  points <- cbind(rep(x, length(y)), y)
  dmvnorm(points, mean = mu, sigma = ht)
}

# 计算X=1处的边缘密度
integrate(myfun_edge_x, lower = -Inf, upper = Inf, x=1, mu=c(0.06,0.03), ht=diag(2))

解析解法（高效准确）

对于联合高斯分布N(μ, Σ)（μ=(μ₁, μ₂)，Σ为2×2协方差矩阵），其边缘密度有解析解：

• X的边缘分布是一维正态分布N(μ₁, Σ₁₁)，其中Σ₁₁是协方差矩阵的(1,1)元素
• Y的边缘分布是一维正态分布N(μ₂, Σ₂₂)，其中Σ₂₂是协方差矩阵的(2,2)元素

直接用dnorm计算即可，无需数值积分，效率更高且结果精确：

# 计算X=1处的边缘密度
dnorm(1, mean = 0.06, sd = sqrt(diag(2)[1]))

# 计算Y=-2处的边缘密度
dnorm(-2, mean = 0.03, sd = sqrt(diag(2)[2]))

你可以对比数值积分和解析解的结果，两者几乎完全一致。

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