optim和glm之間的剩餘標準誤差差異
我嘗試使用擬合甚至R 函數
optim的簡單線性回歸的結果來重現。 參數估計值相同，但殘差估計值和其他參數的標準誤不同，特別是在樣本量較小時。我想這是由於最大似然法和最小二乘法之間計算殘差標準誤差的方式不同（除以 n 或 n-k+1 見下文示例）。 我從網上的閱讀中了解到，優化不是一項簡單的任務，但我想知道是否有可能以簡單的方式重現使用.glm``nls
glm``optim模擬一個小數據集
set.seed(1) n = 4 # very small sample size ! b0 <- 5 b1 <- 2 sigma <- 5 x <- runif(n, 1, 100) y = b0 + b1*x + rnorm(n, 0, sigma)用 optim 估計
negLL <- function(beta, y, x) { b0 <- beta[1] b1 <- beta[2] sigma <- beta[3] yhat <- b0 + b1*x likelihood <- dnorm(y, yhat, sigma) return(-sum(log(likelihood))) } res <- optim(starting.values, negLL, y = y, x = x, hessian=TRUE) estimates <- res$par # Parameters estimates se <- sqrt(diag(solve(res$hessian))) # Standard errors of the estimates cbind(estimates,se) > cbind(estimates,se) estimates se b0 9.016513 5.70999880 b1 1.931119 0.09731153 sigma 4.717216 1.66753138與 glm 和 nls 的比較
> m <- glm(y ~ x) > summary(m)$coefficients Estimate Std. Error t value Pr(>|t|) (Intercept) 9.016113 8.0759837 1.116411 0.380380963 x 1.931130 0.1376334 14.030973 0.005041162 > sqrt(summary(m)$dispersion) # residuals standard error [1] 6.671833 > > summary(nls( y ~ b0 + b1*x, start=list(b0 = 5, b1= 2))) Formula: y ~ b0 + b1 * x Parameters: Estimate Std. Error t value Pr(>|t|) b0 9.0161 8.0760 1.116 0.38038 b1 1.9311 0.1376 14.031 0.00504 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 6.672 on 2 degrees of freedom我可以像這樣重現不同的殘差標準誤差估計：
> # optim / Maximum Likelihood estimate > sqrt(sum(resid(m)^2)/n) [1] 4.717698 > > # Least squares estimate (glm and nls estimates) > k <- 3 # number of parameters > sqrt(sum(resid(m)^2)/(n-k+1)) [1] 6.671833
問題是標準錯誤來自
$$ \hat\sigma^2 (X^\top X)^{-1} $$
在哪裡 $ \hat\sigma^2 $ 是無偏估計量，而不是 MLE。看
summary.lmsummary.lm #R function (object, correlation = FALSE, symbolic.cor = FALSE, #R ...) #R { #R z <- object #R p <- z$rank #R rdf <- z$df.residual #R ... #R Qr <- qr.lm(object) #R ... #R r <- z$residuals #R f <- z$fitted.values #R w <- z$weights #R if (is.null(w)) { #R mss <- if (attr(z$terms, "intercept")) #R sum((f - mean(f))^2) #R else sum(f^2) #R rss <- sum(r^2) #R } #R ... #R resvar <- rss/rdf #R ... #R R <- chol2inv(Qr$qr[p1, p1, drop = FALSE]) #R se <- sqrt(diag(R) * resvar) #R ...這是逆觀察到的 Fisher信息 $ (\beta_0, \beta_1) $ 有條件的 $ \hat\sigma^2 $ . 現在，您計算的反向觀察到的 Fisher 信息是針對三元組的 $ (\beta_0, \beta_1, \sigma) $ . 即，您使用的 MLE $ \sigma $ 而不是無偏估計量。因此，我認為標準誤差應該因因素而異 $ \sqrt{n/(n-3 + 1)} $ 或類似的東西。這是這種情況
set.seed(1) n = 4 # very small sample size ! b0 <- 5 b1 <- 2 sigma <- 5 x <- runif(n, 1, 100) y = b0 + b1*x + rnorm(n, 0, sigma) negLL <- function(beta, y, x) { b0 <- beta[1] b1 <- beta[2] sigma <- beta[3] yhat <- b0 + b1*x return(-sum(dnorm(y, yhat, sigma, log = TRUE))) } res <- optim(c(0, 0, 1), negLL, y = y, x = x, hessian=TRUE) estimates <- res$par # Parameters estimates (se <- sqrt(diag(solve(res$hessian)))) #R [1] 5.690 0.097 1.653 k <- 3 se * sqrt(n / (n-k+1)) #R [1] 8.047 0.137 2.338
為了詳細說明usεr11852的要求，對數似然是
$$ l(\vec{\beta},\sigma) = -\frac{n}{2}\log(2\pi) - n\log{\sigma}
- \frac{1}{2\sigma^2}(\vec{y}-X\vec\beta)^\top(\vec{y}-X\vec\beta) $$
在哪裡 $ X $ 是設計矩陣和 $ n $ 是觀察次數。因此，觀察到的信息矩陣是
$$ -\nabla_{\vec{\beta}}\nabla_{\vec{\beta}}^\top l(\vec{\beta},\sigma) = \frac{1}{\sigma^2}X^\top X $$
現在我們可以插入 MLE 或 $ \sigma $ 如下圖
m <- lm(y ~ x) X <- cbind(1, x) sqrt(sum(resid(m)^2)/n * diag(solve(crossprod(X)))) #R x #R 5.71058285 0.09732149 k <- 3 sqrt(sum(resid(m)^2)/(n-k+1) * diag(solve(crossprod(X)))) #R x #R 8.0759837 0.1376334我們可以對 QR 分解做同樣的
lm事情obj <- qr(X) sqrt(sum(resid(m)^2)/(n-k+1) * diag(chol2inv(obj$qr))) #R [1] 8.0759837 0.1376334所以來回答
我從網上的閱讀中了解到，優化不是一項簡單的任務，但我想知道是否有可能以簡單的方式重現
glm使用optim.那麼您需要在您使用的高斯示例中擴大標準誤差。