如何估計相關觀測值的方差？

假設我們有 n 個觀察值 $ x_i $ （i 從 1 到 n），每個都來自均值為 0 和一些方差分量的正態分佈： $ X_i \sim N(0, \sigma^2) $ . 隨機變量 $ X_i $ s 有一些（假設已知）相關結構。IE： $ Corr(X_i, X_j) = \rho_{ij} $ .

我們應該如何估計 $ \sigma^2 $ ?

相關結構將如何影響估計量，例如 $ \hat \sigma^2 = \frac{\sum(x_i^2)}{n} $ ? 我們能否設計一個更好的估計器（假設相關結構已知）？

讓 $ R = (\rho_{ij}) $ 是相關矩陣，使得協方差矩陣是 $ \Sigma = \sigma^2 R. $ 考慮 $ \mathbf x = (x_1,\ldots, x_n)^\prime $ 成為一個單一的觀察 $ n $ - 均值為零的變量正態分佈和 $ \Sigma $ 協方差。因為對數似然 $ N\ge 1 $ 獨立的這種觀察是，直到一個附加常數（僅取決於 $ N, $ $ n, $ 和 $ R $ ) 由

$$ \Lambda(\sigma) = \sum_{i=1}^N (-n \log \sigma) - \frac{1}{2\sigma^2} \sum_{i = 1}^N\mathbf{x}_i^\prime, R^{-1},\mathbf{x}_i, $$

它有臨界點 $ \sigma\to 0, $ $ \sigma\to\infty, $ 並在任何解決方案

$$ 0 = \frac{\mathrm d}{\mathrm{d}\sigma}\Lambda(\sigma) = -\frac{nN}{\sigma} + \frac{1}{\sigma^3} \sum_{i = 1}^N\mathbf{x}_i^\prime, R^{-1},\mathbf{x}_i. $$

除非表格 $ \sum_{i = 1}^N\mathbf{x}_i^\prime, R^{-1},\mathbf{x}_i $ 為零，有一個唯一的全局最大值

$$ \hat\sigma^2 = \frac{1}{nN} \sum_{i = 1}^N \mathbf{x}_i^\prime, R^{-1},\mathbf{x}_i. $$

這是最大似然估計。 即使在 $ N=1 $ （這是問題中提出的情況）。

直觀地說，這將優於任何忽略假設的相關性的估計 $ R. $ 為了檢查，我們可以計算這個估計值的 Fisher 信息矩陣——但我會把它留給你，部分原因是結果不應該令人信服 $ N $ （最大似然漸近結果可能不適用）。

**為了說明實際發生的情況，**這裡有一千個估計 $ \sigma^2 $ 從實驗 $ n=4 $ 和 $ N=1. $ 的價值 $ \sigma $ 被設置為 $ 1 $ 始終。在這些實驗中，相關性總是

$$ R = \pmatrix{1 &0.3055569 &0.5513377 &0.5100989\ 0.3055569 &1 &0.1240151 &0.09634469\ 0.5513377 &0.12401511 &1 &-0.4209064\ 0.5100989 &0.09634469 &-0.4209064 &1} $$

（一開始隨機生成）。在圖中，最左邊的面板是上述最大似然估計的直方圖；中間面板是使用通常（無偏）方差估計量的估計直方圖；右側面板是兩組估計值的 QQ 圖。斜線是平等線。您可以看到通常的方差估計器往往會產生更多的極值。 它也是有偏差的（由於忽略了相關性）：MLE 的平均值是 0.986 - 驚人地接近真實值 $ \sigma^2 =1^2 =1 $ 而通常估計的平均值僅為 0.791。（我寫“令人驚訝”是因為它是眾所周知的通常的最大似然估計量 $ \sigma^2, $ 在不涉及相關性的情況下，具有順序偏差 $ 1/(nN), $ 在這種情況下相當大。）

您可以通過修改、、、、和隨機數種子R的值來試驗生成這些數字的代碼。n``sigma``N``n.sim``Rho``17

f <- function(x, Rho) { # The MLE of sigma^2 given data `x` and correlation `Rho`
 S <- solve(Rho)
 sum(apply(x, 1, function(x) x %*% S %*% x)) / length(c(x))
}
n <- 4
sigma <- 1
N <- 1
n.sim <- 1e3
set.seed(17)
#
# Generate a random correlation matrix. Larger values of `d` yield more
# spherical matrices in general.
#
d <- 1
Rho <- cor(matrix(rnorm(n*(n+d)), ncol=n))
(ev <- eigen(Rho, only.values=TRUE)$values) 
#
# Run the experiments.
#
library(MASS)
sim <- replicate(n.sim, {
 x <- matrix(mvrnorm(N, rep(0,n), sigma^2 * Rho), N)
 c(f(x, Rho),  var(c(x)))
})
(rowMeans(sim))
#
# Plot the results.
#
par(mfrow=c(1,3))
hist(sim[1,], col=gray(.93), xlab="Estimate", 
    main=expression(paste("Histogram of Estimates of ", sigma^2)))
abline(v = sigma^2, col="Red", lwd=2)

hist(sim[2,],col=gray(.93), xlab="Estimate",
    main=expression(paste("Histogram of Independent Estimates of ", sigma^2)))
abline(v = sigma^2, col="Red", lwd=2)

y1 <- sort(sim[1,])
y2 <- sort(sim[2,])
plot(y1, y2,asp=1, xlab="Correlation-based estimate", ylab="Independent estimate")
abline(0:1, col="Red", lwd=2)
par(mfrow=c(1,1))

引用自：https://stats.stackexchange.com/questions/535727