Первоначально я хотел бы сделать два предложения, которые могли бы сделать проблему проще для кода. Во-первых, пересмотреть схему данных, чтобы каждый год представлял собой уникальную строку, и каждая группа является уникальным столбцом. Во-вторых, поскольку когорты обрабатываются математически независимо друг от друга, держите их отдельно на время, по крайней мере, до тех пор, пока не будет построено ядро кода. Помещенный цикл вокруг этого позже, который циклически проходит через них. В первом блоке кода есть две матрицы, одна с наблюдаемыми данными, и одна, которая будет собирать предсказанные данные.
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 воссоздает память с помощью функции замены).
Спасибо, WIll. Это именно тот механизм, который я искал. Моя ошибка заключалась в том, чтобы думать о широких формах данных при кодировании модельных уравнений в петли. Широкоформатное преобразование немного привыкает, но в конце концов должно заплатить. – andrey