2013-03-12 4 views
1

Работа в 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 и представляют соответствующий размер каждой группы (А, В, С) в момент времени Т = 00. Они добавляют до 1 по строке (A00 + B00 + C00 = 1). Параметры Вкладка и Тас используются для прогнозирования распространенности в момент времени Т + 1, используя некоторую математическую модель, например

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

Функция для вычисления предсказанных значений в момент времени Т + 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, например). Заранее спасибо!

ответ

2

Первоначально я хотел бы сделать два предложения, которые могли бы сделать проблему проще для кода. Во-первых, пересмотреть схему данных, чтобы каждый год представлял собой уникальную строку, и каждая группа является уникальным столбцом. Во-вторых, поскольку когорты обрабатываются математически независимо друг от друга, держите их отдельно на время, по крайней мере, до тех пор, пока не будет построено ядро ​​кода. Помещенный цикл вокруг этого позже, который циклически проходит через них. В первом блоке кода есть две матрицы, одна с наблюдаемыми данными, и одна, которая будет собирать предсказанные данные.

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. Я использую 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 loop немного чище, чем когда каждая пара элементов была явно умножена и суммирована. Но мы не должны рассматривать каждую пару строк/столбцов отдельно. Все три столбца ob матрицы могут работать на всех трех столбцов tt матрицы одновременно:

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

Это должно быть намного быстрее, чем даже в предыдущей версии, так как внутренняя память системы R является не воссоздает матрицу три раза для каждой строки - только один раз в строке. Чтобы уменьшить это до одного раза на матрицу, используйте функцию apply, а затем транспонируйте матрицу, если это соответствует вашей цели. Наконец, обратите внимание, что строки представляют разные годы, чем pred (т. Е. Строка i-1 здесь такая же, как строка i в pred).

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

Для размещения когорт, возможно, вы могли бы объявить список с тремя элементами (для 1980, 1981 и 1982 когорты).Каждый элемент будет уникальной матрицей ob. И создайте второй список для уникальной матрицы pred. Или, может быть, использовать трехмерные матрицы (но это может быть больше налогов, когда R воссоздает память с помощью функции замены).

+0

Спасибо, WIll. Это именно тот механизм, который я искал. Моя ошибка заключалась в том, чтобы думать о широких формах данных при кодировании модельных уравнений в петли. Широкоформатное преобразование немного привыкает, но в конце концов должно заплатить. – andrey

Смежные вопросы