1

我在 Rstudio 中使用集成卡尔曼滤波器来查看 Sir 模型。我有 i = 100 个整体成员的预测值 (xf) ^ i 和 I 的测量值 (z)。我现在想通过以下方式组合它们：

  • 对于时间步 1 到 10，我想使用预测值。
  • 然后在 t = 10 时，我做一个分析值： ((xa) ^ i = (xf) ^ i + k (z - M * (xf) ^ i + v ^ i)，其中 k 是权重矩阵， v ^ i = 随机测量误差，我在 t = 10 时输入 (xa) ^ i。
  • 现在我想用这些 t=10 的值作为新的初始值来计算 t = 11 tm t = 20 的预测值。对于 t = 20 和 t = 10 会发生同样的事情。

但是，以相同的方式计算预测值，但使用新的初始值似乎不起作用。我的代码如下。有人知道怎么做这个吗？谢谢！

#kalman gain vaste waarden 
#matrix M, 1 bij variabele die we bekijken
#aangeven welke digit 1 moet worden M[,i]=1, S: i=1, I: i=2, R: i=3
M <- matrix(data=0, nrow=1, ncol=3)
M[,2] = 1
zvec <- as.vector(read_excel("z.xlsx", range = "B1:B101")) 
z <- t(zvec)
Rv = z/20 #is variantie van de vi trekking, var = sd^2
vi <- rnorm(n=1001, mean=0, sd = Rv) 
tijdstap <- seq(0, 999, stap) #aangeven om de hoeveel t's er een meting
xfmat <- matrix(0, 3, 100) #voor losse t xfS, xfI en xfR samenvoegen tot forecast matrix
xkl <- list()
xa <- list()

xfS <- matrix(0, 1000, 100)
xfI <- matrix(0, 1000, 100)
xfR <- matrix(0, 1000, 100)
for(u in 1:10){
  xfS[u,] = SS[u,]
  xfI[u,] = II[u,]
  xfR[u,] = RR[u,]
}

for(t in 10:999){
  if(t %in% tijdstap){
    xfmat[1,] = xfS[t,]
    xfmat[2,] = xfI[t,]
    xfmat[3,] = xfR[t,]
    
    k = (covarmatrix[[t]] %*% t(M)) %*% (1/(M %*% covarmatrix[[t]] %*% t(M) + Rv[t/10]))
    xa[[t]] = xfmat + k%*%(z - M%*%xfmat + vi[t/10]) 
    xfS[t,] = xa[[t]][1,] 
    xfI[t,] = xa[[t]][2,]
    xfR[t,] = xa[[t]][3,]
    
    #using xa as new initial values for xf --> this doesn't seem to work
    for(i in 1:100){
      xkl[[i]] <- xkfun(xa[[t]][1,i], xa[[t]][2,i], xa[[t]][3,i], rnorm(1, mean = 1, sd = 0.5), 0.1, )
      xfS[(t+1):(t+10),i] = xkl[[i]][,1]
      xfI[(t+1):(t+10),i] = xkl[[i]][,2]
      xfR[(t+1):(t+10),i] = xkl[[i]][,3]
    }
    
    #new mean and covariance matrix for new values
    gemnew <- gemfun(xfS[(1:(t+10)),], xfR[(1:(t+10)),], xfR[(1:(t+10)),])
    gemnew <- gemfun(xfS, xfI, xfR)
    covarmatrix <- covarmatfun(xfS, xfI, xfR, gemnew, (t+10))
  }
}
4

0 回答 0