我正在做一個研究專案，不久前我在數學堆疊交換上問了這個問題，在那里我正在尋找一種方法來計算給定轉換矩陣的收入期間變化的方差，其中每個狀態對應到給定向量中收入的對數水平。我想計算一個人在每個州開始的收入變化在 n 個時期內的方差是多少。我的狀態空間由 11 個狀態組成，所以我希望最終得到一個包含 11 個不同方差的向量。當我問這個問題時，我得到了一個滿意的答案，但是當我嘗試在 RI 中撰寫代碼時遇到了一些問題，希望能得到幫助。
我創建了這段代碼來計算方差：
install.packages("expm")
library(expm)
# creating standard basis vectors
e <- function(i) {
e_i = rep(0, length(alpha))
e_i[i] = 1
return(e_i)
}
# compute variances
p2p_variance = function(n, alpha, P) {
variance = list()
pi_n = list()
for (i in 1:length(alpha)) {
pi_n[[i]] = e(i) %*% (P %^% n)
beta = (t(alpha) - t(alpha)[i])^2
variance[[i]] = (pi_n[[i]] %*% t(beta)) - (((pi_n[[i]] %*% alpha) - alpha[i]) %^% 2)
}
return(t(variance))
}
對于我的 alpha（收入對數水平向量）和 P（轉換矩陣）的值，我使用：
alpha = c(3.4965, 3.5835, 3.6636, 3.7377, 3.8067, 3.8712, 3.9318, 3.9890, 4.0431, 4.0943, 4.1431)
P = rbind(c(0.9004, 0.0734, 0.0203, 0.0043, 0.0010, 0.0003, 0.0001, 0.0001, 0.0000, 0.0000, 0.0000),
c(0.3359, 0.3498, 0.2401, 0.0589, 0.0115, 0.0026, 0.0007, 0.0003, 0.0001, 0.0001, 0.0000),
c(0.1583, 0.1538, 0.3931, 0.2346, 0.0481, 0.0090, 0.0021, 0.0007, 0.0003, 0.0001, 0.0001),
c(0.0746, 0.0609, 0.1600, 0.4368, 0.2178, 0.0397, 0.0073, 0.0019, 0.0006, 0.0002, 0.0001),
c(0.0349, 0.0271, 0.0559, 0.1724, 0.4628, 0.2031, 0.0344, 0.0067, 0.0018, 0.0006, 0.0003),
c(0.0155, 0.0122, 0.0230, 0.0537, 0.1817, 0.4870, 0.1860, 0.0316, 0.0066, 0.0018, 0.0009),
c(0.0066, 0.0054, 0.0100, 0.0204, 0.0529, 0.1956, 0.4925, 0.1772, 0.0307, 0.0064, 0.0023),
c(0.0025, 0.0023, 0.0043, 0.0084, 0.0186, 0.0530, 0.2025, 0.4980, 0.1760, 0.0275, 0.0067),
c(0.0009, 0.0009, 0.0017, 0.0035, 0.0072, 0.0168, 0.0490, 0.2025, 0.5194, 0.1721, 0.0260),
c(0.0003, 0.0003, 0.0007, 0.0013, 0.0029, 0.0061, 0.0142, 0.0430, 0.2023, 0.5485, 0.1804),
c(0.0001, 0.0001, 0.0002, 0.0003, 0.0008, 0.0017, 0.0032, 0.0068, 0.0212, 0.1079, 0.8578))
例如，呼叫p2p_variance(100, alpha, P)（計算 100 個周期的方差）會產生以下方差向量：
0.04393012 0.04091066 0.03856503 0.03636202 0.03472286 0.03331921 0.03213084 0.03068901 0.03143765 0.03255994 0.03522346
這似乎是合理的。但是，如果我運行p2p_variance(1000, alpha, P)，它會導致：
0.06126449 0.03445073 0.009621497 -0.01447615 -0.03652425 -0.05752316 -0.07753646 -0.09726683 -0.1134972 -0.1287498 -0.141676
這顯然是不正確的，因為我們不能有負方差。我無法弄清楚為什么簡單地將 n 增加到 1000 會導致這里出現負方差。我很可能錯誤地編碼了我的 p2p_variance 函式，但我終生無法找到問題。或者我用來發現這些差異的程序是否存在某種缺陷？如果有人可以查看此代碼并幫助我診斷問題，我將不勝感激
uj5u.com熱心網友回復：
您的方差函式回傳差異，如果您想要絕對值（方差），只需將其包裝在里面，abs()如下所示：
p2p_variance = function(n, alpha, P) {
variance = list()
pi_n = list()
for (i in 1:length(alpha)) {
pi_n[[i]] = e(i) %*% (P %^% n)
beta = (t(alpha) - t(alpha)[i])^2
variance[[i]] = abs((pi_n[[i]] %*% t(beta)) - (((pi_n[[i]] %*% alpha) - alpha[i]) %^% 2))
}
return(t(variance))
}
p2p_variance(1000, alpha, P)
輸出：
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
[1,] 0.06126449 0.03445073 0.009621497 0.01447615 0.03652425 0.05752316 0.07753646 0.09726683 0.1134972 0.1287498 0.141676
轉載請註明出處，本文鏈接：https://www.uj5u.com/yidong/446276.html