时间序列预测函数

Question

在 R 中工作。我想使用初始值和一组转换参数来预测患病率的时间序列。对于如下结构的数据

 cohort <- c(1980,1981,1982)
 A00 <- c(.15, .2,.4)
 B00 <- c(.25, .3, .4) 
 C00 <-c(.6, .5,.2)
 Tab<-c(.6,.5,.4)
 Tac<-c(.2,.25,.35)
 ds <- data.frame(cohort,A00,B00,C00,Tab,Tac)
 print (ds)

  cohort  A00  B00 C00 Tab  Tac
1   1980 0.15 0.25 0.6 0.6 0.20
2   1981 0.20 0.30 0.5 0.5 0.25
3   1982 0.40 0.40 0.2 0.4 0.35

A00、B00 和 C00 列中的初始值表示每个组 (A,B,C) 在时间 t=00 的相关大小。它们在整行中加起来为 1 (A00+B00+C00=1)。参数 Tab 和 Tac 用于使用一些数学模型预测时间 t+1 的患病率，例如

A01   = df$A00 -df$Tab +df$Tac.

计算t+1时刻预测值的函数是

 forecast<- function( df ) {
  dsResult <- data.frame(
    cohort= df$cohort,
    A01   = df$A00 -df$Tab +df$Tac ,    
    B01   = df$B00 -df$Tab +df$Tac,    
    C01  =  df$C00 -df$Tab +df$Tac    

  )
  dsResult<- merge(df,dsResult,by="cohort")
  return( dsResult)
}
new<-forecast(ds)

并产生以下结果

  cohort  A00  B00 C00 Tab  Tac   A01   B01  C01
1   1980 0.15 0.25 0.6 0.6 0.20 -0.25 -0.15 0.20
2   1981 0.20 0.30 0.5 0.5 0.25 -0.05  0.05 0.25
3   1982 0.40 0.40 0.2 0.4 0.35  0.35  0.35 0.15

非常感谢您帮助我学习如何编写一个循环来循环预测所需的年数（例如，对于 1:7 中的 t）。提前致谢！

Answer 1

最初，我想提出两个建议，以使问题更易于编码。首先，修改数据模式，使每一年都是唯一的行，每个组是唯一的列。其次，由于队列在数学上彼此独立，因此暂时将它们分开，至少在构建代码的内核之前。稍后围绕这个循环，循环遍历它们。在第一个代码块中，有两个矩阵，一个用于观察数据，一个用于收集预测数据。

yearCount <- 7 #Declare the number of time points.
groupCount <- 3 #Declare the number of groups.

#Create fake data that sum to 1 across rows/times.
ob <- matrix(runif(yearCount*groupCount), ncol=groupCount)
ob <- ob / apply(ob, 1, function( x ){ return( sum(x) )})

#Establish a container to old the predicted values.
pred <- matrix(NA_real_, ncol=groupCount, nrow=yearCount)

t12<-.5; t13<-.2; t11<-1-t12-t13 #Transition parameters from group 1
t21<-.2; t23<-.4; t22<-1-t21-t23 #Transition parameters from group 2
t31<-.3; t32<-.1; t33<-1-t31-t32 #Transition parameters from group 3

for( i in 2:yearCount ) {
  pred[i, 1] <- ob[i-1, 1]*t11 + ob[i-1, 2]*t21 + ob[i-1, 3]*t31
  pred[i, 2] <- ob[i-1, 1]*t12 + ob[i-1, 2]*t22 + ob[i-1, 3]*t32
  pred[i, 3] <- ob[i-1, 1]*t13 + ob[i-1, 2]*t23 + ob[i-1, 3]*t33
}

#Calculate the squared errors
ss <- (pred[-1, ] - ob[-1, ])^2 #Ignore the first year of data

在循环内部，您可能会注意到矩阵乘法的熟悉结构。每行可以使用内积稍微压缩（即，将ob 矩阵的一行相乘，然后与ts 的一个“列”相加。我使用t12 与@ 略有不同987654325@ 在您的帖子中；这是在给定时间点从第 1 组转换到第 2 组的概率。

#Create transition parameters that sum to 1 across rows/groups.
tt <-  matrix(runif(groupCount*groupCount), ncol=groupCount)
tt <- tt / apply(tt, 1, function( x ){ return( sum(x) )})

假设tt 矩阵是之前定义的，而不是t11,...,t33 的单独变量。

for( i in 2:yearCount ) {
  pred[i, 1] <- ob[i-1, ] %*% tt[, 1] 
  pred[i, 2] <- ob[i-1, ] %*% tt[, 2]
  pred[i, 3] <- ob[i-1, ] %*% tt[, 3]
}

循环的内容比每个元素对显式相乘和相加时稍微干净一些。但是我们不必单独处理每一行/列对。 ob矩阵的所有三列都可以被tt矩阵的所有三列同时操作：

for( i in 2:yearCount ) {
  pred[i, ] <- ob[i-1, ] %*% tt
}

这应该比以前的版本快得多，因为 R 的内部存储系统不会为每行重新创建矩阵三次 - 每行仅一次。要将其减少到每个矩阵一次，请使用apply 函数，然后在适合您的目的时转置矩阵。最后，请注意，这些行表示的年份与pred 不同（即此处的第 i-1 行与pred 中的第 i 行相同）。

predictionWIthExtraYear <- t(apply(ob, 1, FUN=function(row){row %*% tt}))

为了容纳同类群组，也许您可​​以声明一个包含三个元素的列表（针对 1980、1981 和 1982 年的同类群组）。每个元素都是一个唯一的ob 矩阵。并为唯一的pred 矩阵创建第二个列表。或者可以使用三维矩阵（但是当 R 使用替换函数重新创建内存时，这可能会更加繁重）。