QR 분해를 배우려고하지만 전통적인 행렬 계산을 사용하지 않고 beta_hat의 분산을 얻는 방법을 알 수 없습니다. 나는 iris 데이터 세트와 함께 연습하고, 여기에 내가 지금까지 가지고 무엇을 해요 : 당신의 도움에 대한R에서 QR 분해를 사용하여 최소 제곱 추정량의 분산을 계산하는 방법?
y<-(iris$Sepal.Length)
x<-(iris$Sepal.Width)
X<-cbind(1,x)
n<-nrow(X)
p<-ncol(X)
qr.X<-qr(X)
b<-(t(qr.Q(qr.X)) %*% y)[1:p]
R<-qr.R(qr.X)
beta<-as.vector(backsolve(R,b))
res<-as.vector(y-X %*% beta)
감사합니다!
잔류 정도 n - p이 (더 정확하게 n - qr.X$rank. 이하에서, I는 qr.X$rank 대신 p를 사용한다.) 때문에 추정 분산 따라서
se2 <- sum(res^2)/(n - qr.X$rank)
이다 이며 beta_hat의 분산 공분산 행렬은 단지
chol2inv(R) * se2
# [,1] [,2]
#[1,] 0.22934170 -0.07352916
#[2,] -0.07352916 0.02405009
매우 자주, 우리는 단지 변이 분산, 즉이 완전한 분산 - 공분산 행렬의 대각선만을 원한다. 이 경우 chol2inv을 먼저 호출하여 전체 공분산 행렬을 명시 적으로 형성 할 필요가 없습니다. 우리는 훨씬 낮은 계산 비용으로 직접 대각선을 계산할 수 있습니다 :
Rinv <- backsolve(R, diag(qr.X$rank))
rowSums(Rinv^2) * se2
# [1] 0.22934170 0.02405009
이의이 lm와 비교하여 정확성을 확인하자
fit <- lm(Sepal.Length ~ Sepal.Width, iris)
vcov(fit)
# (Intercept) Sepal.Width
#(Intercept) 0.22934170 -0.07352916
#Sepal.Width -0.07352916 0.02405009
예, 모든 올바른!
마지막으로, 게시 코드에 대한 몇 가지 코멘트 :
대신
는b <- (t(qr.Q(qr.X)) %*% y)[1:p]
당신은 당신이없는
b <- qr.qty(qr.X, y)[1:qr.X$rank]
을 사용할 수 있습니다를 추출하십시오. backsolve에 대한; qr.X$qr을 사용하는 것만으로는 충분하다 :
beta <- as.vector(backsolve(qr.X$qr,b))
Rinv <- backsolve(qr.X$qr, diag(qr.X$rank))
가 여기에 완전한 세션 : 대답을 수락
## data
y <- iris$Sepal.Length
x <- iris$Sepal.Width
## add intercept to model matrix
X <- cbind(1,x)
## QR factorization of model matrix
qr.X <- qr(X)
## get estimated coefficient
b <- qr.qty(qr.X, y)
beta <- as.vector(backsolve(qr.X$qr, b))
## compute residuals (I don't use `qr.resid` and `qr.fitted`, though I can)
res <- as.vector(y - X %*% beta)
## residual standard error
se2 <- sum(res^2)/(nrow(X) - qr.X$rank)
## full variance-covariance matrix
chol2inv(qr.X$qr) * se2
@ZheyuanLi입니다 대부분의 답변이 작동하는지 표시 및을 제공에 대해 조금 덜에 대한 응답자가 투자 한 노력에 대한 보상. 반드시 전문적인 행동에 관한 것은 아닙니다. 질문을 (완전히) 해결하지 못했을 때 대답을 수락하지 않으면 전문적인 행동으로 분류 될 수 있습니다. – Jaap
@ ZheyuanLi 안녕하세요, 저는 길 위에 있었고 모바일에서 수락 할 수 없었습니다 .... 감사합니다! – Nudibranch