Р: Может ли матрица заполняться в цикле путем индексации строк и столбцов?
Я биолог и новичок в R, и я учусь создавать простую популяционную модель.
Итак, у меня есть популяционная матрица ("поп") из 30 возрастных классов женщин (1:4 - не заводчики, 5:30 - заводчики), которые будут моделироваться в течение 100 лет.
pop <- matrix(0,30,100)
Затем я заполняю эту матрицу 3 молодыми взрослыми женщинами.
pop[5, 1] <- 3
Затем я хочу выполнить это в течение 100 лет со стохастичностью, чтобы увидеть, как это население делает со временем. (Я не заполнил их, но они вам не нужны, все они имеют разные вероятности выживания по сравнению с половозрелыми взрослыми.)
for (t in 1:100) { # Edited y to t, typing error!
pop[1,t+1] <- rbinom(1,colSums(pop[5:30, t]), b/2)
pop[5, t+1] <- rbinom(1, pop[4, t], s2)
pop[6, t+1] <- rbinom(1, pop[5, t], s2)
.....
pop[30, t+1] <- rbinom(1, pop[29, t], s2)
}
Итак, мой вопрос: есть ли способ заполнить эту матрицу без необходимости явно писать 30 строк кода? Поскольку строки 5 - 30 все будут одинаковыми, но даже после 5 часов (буквально) веб-поиска и чтения R вручную, я не могу найти способ индексировать строки, что, по-видимому, и здесь необходимо.
Любое понимание здесь приветствуется, в том числе другой способ моделирования этой популяции.
1 ответ
Поскольку rbinom
векторизация, вы должны быть в состоянии заменить
pop[1,t+1] <- rbinom(1,colSums(pop[5:30, t]), b/2)
pop[5, t+1] <- rbinom(1, pop[4, t], s2)
pop[6, t+1] <- rbinom(1, pop[5, t], s2)
.....
pop[30, t+1] <- rbinom(1, pop[29, t], s2)
с
pop[,t+1] <- rbinom(27,size=c(colSums(pop[5:30,t]),pop[4:29,t]),
prob=c(b/2,rep(s2,26)))
или что-то в этом роде (26 и 27 могут быть не правы, я, возможно, не учел)
Тем не менее, я действительно думаю, что вы должны сделать что-то вроде этого:
pop[1,t+1] <- rpois(1,b/2*sum(pop[5:30,t]))
pop[-1,t+1] <- rbinom(29,size=pop[1:29,t]),
prob=c(rep(s1,4),rep(s2,26)))
Это предполагает, что люди в последнем возрастном классе падают с края света и умирают... Я использую rpois
здесь, чтобы разрешить воспроизведение Пуассона, а не бином, хотя вы могли бы использовать rbinom(1,size=sum(pop[5:30,t])/2,prob=p)
если вы действительно хотите настаивать на том, чтобы у самок было максимум 1 потомство (с вероятностью p)