在R中手动计算OLS估计器遇维度问题及自动化实现求助

一、维度错误的原因及修正 #

你的代码存在两个核心问题导致维度报错：

1. 矩阵乘法运算符误用：R中*是元素逐乘，矩阵叉乘必须用%*%。原代码t(x)*x生成的是与x同维度的矩阵，而OLS公式需要的是t(x)%*%x（1000×1000的方阵），这直接导致后续求逆和乘法的维度不匹配。
2. 求逆函数的选择：建议用R基础包的solve()做矩阵求逆；若需处理奇异矩阵，可改用MASS::ginv()（广义逆）。

修正后的单变量计算代码：

library(MASS) # 使用广义逆时加载

nr = 100
nc = 1000 
x <- matrix(rnorm(nr * nc, mean = 1, sd = 1), nrow = nr)
epsilon <- matrix(rnorm(nr * nc, mean = 0, sd = 1), nrow = nr)
k <- c(1,2,4,8)

# 用scale简化中心化标准化（与手动公式等价）
eta1 <- scale(epsilon^1, center = TRUE, scale = TRUE)
eta2 <- scale(epsilon^2, center = TRUE, scale = TRUE)
eta4 <- scale(epsilon^4, center = TRUE, scale = TRUE)
eta8 <- scale(epsilon^8, center = TRUE, scale = TRUE)

y1 <- x + eta1
y2 <- x + eta2
y4 <- x + eta4
y8 <- x + eta8

# 修正矩阵乘法，用%*%
xtx <- t(x) %*% x
beta1 <- solve(xtx) %*% t(x) %*% y1
# 奇异矩阵时替换为：beta1 <- ginv(xtx) %*% t(x) %*% y1

二、基于k值的自动化批量计算 #

不用手动生成每个eta和beta，用lapply遍历k向量实现自动化，同时预计算t(x)%*%x提升效率：

set.seed(123) # 设置随机种子，保证结果可复现
nr = 100
nc = 1000 
x <- matrix(rnorm(nr * nc, mean = 1, sd = 1), nrow = nr)
epsilon <- matrix(rnorm(nr * nc, mean = 0, sd = 1), nrow = nr)
k_vec <- c(1,2,4,8)

# 预计算xtx，避免重复计算
xtx <- t(x) %*% x

# 遍历每个k值，生成对应的beta
beta_results <- lapply(k_vec, function(current_k) {
  # 计算标准化后的eta
  eta <- scale(epsilon^current_k, center = TRUE, scale = TRUE)
  # 构造y矩阵
  y <- x + eta
  # 计算OLS估计量
  beta <- solve(xtx) %*% t(x) %*% y
  # 返回k值和对应的beta矩阵
  list(k = current_k, beta_matrix = beta)
})

# 查看k=1对应的beta矩阵
beta_results[[1]]$beta_matrix

# 可选：将所有beta合并为一个矩阵，列名对应k值
beta_combined <- do.call(cbind, lapply(beta_results, function(res) res$beta_matrix))
colnames(beta_combined) <- paste0("k=", k_vec)

关键说明 #

• scale()函数完全等价于你手动编写的中心化标准化公式，代码更简洁且不易出错。
• 预计算xtx可大幅减少重复计算量，因为所有k值对应的xtx是相同的。
• 若xtx存在线性相关性（奇异矩阵），将solve()替换为MASS::ginv()即可解决报错。

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