我试图在R中计算P任意N x J矩阵的投影矩阵S:
P = S (S'S) ^ -1 S'
我一直试图用以下功能执行此操作:
P <- function(S){
output <- S %*% solve(t(S) %*% S) %*% t(S)
return(output)
}
但是当我使用它时,我得到的错误看起来像这样:
# Error in solve.default(t(S) %*% S, t(S), tol = 1e-07) :
# system is computationally singular: reciprocal condition number = 2.26005e-28
我认为这是数值下溢和/或不稳定的结果在许多地方一样讨论R-帮助和这里,但我没有足够的经验使用SVD或QR分解来解决这个问题,要不然就把这个现有的代码到行动.我也尝试了建议的代码,即将solve写成一个系统:
output <- S %*% solve (t(S) %*% S, t(S), tol=1e-7)
但它仍然无效.任何建议,将不胜感激.
我很确定我的矩阵应该是可逆的并且没有任何共线性,只是因为我尝试用正交虚拟变量矩阵进行测试,但它仍然不起作用.
另外,我想将它应用于相当大的矩阵,所以我正在寻找一个简洁的通用解决方案.
我正在对具有多个预测变量的线性模型的预测值求和，如下面的示例所示，并希望计算该总和的组合方差，标准误差和可能的置信区间。
lm.tree <- lm(Volume ~ poly(Girth,2), data = trees)
假设我有一组Girths：
newdat <- list(Girth = c(10,12,14,16)
为此，我想预测总数Volume：
pr <- predict(lm.tree, newdat, se.fit = TRUE)
total <- sum(pr$fit)
# [1] 111.512
如何获得方差total？
我正在尝试扩展lwr()包的功能McSptial，它适合将加权回归作为非参数估计。在函数的核心中lwr()，它使用以下方法反转矩阵solve()，它使用QR 分解而不是 QR 分解来我想更改它，但无法弄清楚如何从 QR 分解中获取帽子矩阵（或其他导数）。
有数据：
set.seed(0); xmat <- matrix(rnorm(500), nrow=50) ## model matrix
y <- rowSums(rep(2:11,each=50)*xmat) ## arbitrary values to let `lm.wfit` work
w <- runif(50, 1, 2) ## weights
该lwr()功能如下：
xmat2 <- w * xmat
xx <- solve(crossprod(xmat, xmat2))
xmat1 <- tcrossprod(xx, xmat2)
vmat <- tcrossprod(xmat1)
我需要的值，例如：
sum((xmat[1,] %*% xmat1)^2)
sqrt(diag(vmat))
目前我使用reg <- lm.wfit(x=xmat, y=y, w=w)但无法设法恢复在我看来是帽子矩阵（xmat1）的东西。
rich.main3是R中的线性模型。我了解列表的其余元素，但我不知道这qraux是什么。文档指出它是
一个长度为ncol（x）的向量，其中包含有关\ bold {Q}“的附加信息。
这意味着什么附加信息？
str(rich.main3$qr)
qr : num [1:164, 1:147] -12.8062 0.0781 0.0781 0.0781 0.0781 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:164] "1" "2" "3" "4" ...
.. ..$ : chr [1:147] "(Intercept)" "S2" "S3" "x1" ...
..- attr(*, "assign")= int [1:147] 0 1 1 2 3 4 5 6 7 8 ...
..- attr(*, "contrasts")=List of 3
.. ..$ S : chr "contr.treatment"
.. ..$ ID : chr "contr.treatment"
.. ..$ …我使用glmnet R包运行Ridge回归.我注意到我从glmnet::glmnet函数中获得的系数与通过定义计算系数得到的系数不同(使用相同的lambda值).有人可以解释一下为什么吗？
数据(两者:响应Y和设计矩阵X)都是按比例缩放的.
library(MASS)
library(glmnet)
# Data dimensions
p.tmp <- 100
n.tmp <- 100
# Data objects
set.seed(1)
X <- scale(mvrnorm(n.tmp, mu = rep(0, p.tmp), Sigma = diag(p.tmp)))
beta <- rep(0, p.tmp)
beta[sample(1:p.tmp, 10, replace = FALSE)] <- 10
Y.true <- X %*% beta
Y <- scale(Y.true + matrix(rnorm(n.tmp))) # Y.true + Gaussian noise
# Run glmnet
ridge.fit.cv <- cv.glmnet(X, Y, alpha = 0)
ridge.fit.lambda <- ridge.fit.cv$lambda.1se
# Extract coefficient values for …我只是在拼命尝试让乔列斯基分解发挥作用，以模拟相关的价格变动。
我使用以下代码：
cormat <- as.matrix(read.csv("http://pastebin.com/raw/qGbkfiyA"))
cormat <- cormat[,2:ncol(cormat)]
rownames(cormat) <- colnames(cormat)
cormat <- apply(cormat,c(1,2),FUN = function(x) as.numeric(x))
chol(cormat)
#Error in chol.default(cormat) :
# the leading minor of order 8 is not positive definite
cholmat <- chol(cormat, pivot=TRUE)
#Warning message:
# In chol.default(cormat, pivot = TRUE) :
# the matrix is either rank-deficient or indefinite
rands <- array(rnorm(ncol(cholmat)), dim = c(10000,ncol(cholmat)))
V <- t(t(cholmat) %*% t(rands))
#Check for similarity
cor(V) - cormat ## Not all zeros!
#Check the standard deviations
apply(V,2,sd) …我有一个要适合数据的公式列表，而不是运行一个循环，出于性能考虑，我想立即执行此操作。估算应该仍然是分开的，我不是要估算SUR或其他任何值。下面的代码做我想要的
x <- matrix(rnorm(300),ncol=3)
y <- x %*% c(1,2,3)+rnorm(100)
formulae <-list(y~x[,1],
y~x[,2],
y~x[,1] + x[,2])
lapply(formulae,lm)
不幸的是，formulae随着增加长度的增加，这变得有些慢了，有没有办法真正将其向量化？
如果有帮助，lm我唯一关心的结果就是系数和一些标准误差。