按時間序列預測的函數

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和在時間t = 00表示各組（A，B，C）的相關尺寸。它們在整行上加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， 例如）。提前致謝！

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矩陣的一行相乘，然後與t s的一個「列」求和。我使用t12與您的文章中的Tab略有不同;這是在給定時間點從第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] 
} 

Th e循環的內容比每個元素對被明確相乘和求和時稍微乾淨一些。但是我們不必分別處理每個行/列對。在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用替換函數重新創建內存時，這可能會導致更多的稅收）。