关于R中Dagenais-Dagenais高阶矩估计z计算代码的正确性验证

你的代码中part2部分存在错误，不符合公式要求，问题出在矩阵均值计算和后续运算逻辑上，具体分析和修正方案如下：

公式拆解与原代码问题 #

原公式里的E[(x'x)/T] * I[K]，是K×K二阶矩矩阵与K阶单位矩阵的Hadamard逐元素乘积，结果是仅保留二阶矩矩阵对角线元素的对角矩阵。而你代码里的mean((t(x.mat) %*% x.mat) / nrow(x.mat) * diag(ncol(x.mat)))会把这个对角矩阵的所有元素（包括大量0）求平均值，得到一个标量后再与x.mat做矩阵乘法，完全偏离了公式意图。

修正后的代码 #

# y.vec 30 by 1
# x.mat 30 by 3
set.seed(123)
Y <- rt(n = 30, df = 4)
X <- matrix(data = rt(n = 90, df = 4),
            nrow = 30, ncol = 3,
            dimnames = list(NULL, c("X1", "X2", "X3")))
# Y and X in mean deviations
y.vec <- Y - mean(Y)  # 简化中心化逻辑，无需转data.frame调用apply
x.mat <- apply(X, 2, function(a) a - mean(a))

# 计算z
part1 <- x.mat * x.mat * x.mat  # 正确的Hadamard三次方乘积

# 修正part2的计算
T <- nrow(x.mat)
K <- ncol(x.mat)
# 计算样本二阶矩矩阵
S <- crossprod(x.mat) / T
# 生成S与单位矩阵的Hadamard乘积（仅保留对角线的对角矩阵）
D <- diag(diag(S))
part2 <- 3 * x.mat %*% D

z <- part1 - part2

额外优化说明 #

• 中心化y.vec可直接写成Y - mean(Y)，无需转成data.frame调用apply，逻辑更简洁。
• 使用crossprod(x.mat)替代t(x.mat) %*% x.mat，R内部实现更高效。
• 等价简洁写法：由于x.mat %*% D等价于x.mat与S的对角线元素逐元素相乘，part2也可写成：

  part2 <- 3 * x.mat * rep(diag(S), each = T)
  

  利用R的广播特性，结果完全一致。

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