R中的长时间序列函数的简化/浓缩

我有一个很长的时间序列循环，我想简化/压缩。我试图用随机二项式分布模拟一个牛群在十年期间（每月间隔）的产犊。该功能始于假设牛已被公牛覆盖。每个变量都受到之前的影响。变量如下：R中的长时间序列函数的简化/浓缩

G1：G9妊娠每个月。 MC1：MC7与小牛7个月的母亲，然后小牛断奶后。 Rest1：Rest6休息期再次被公牛覆盖。 DeadCows基于死亡率。基于受孕率的NPreg非怀孕母牛。

输入： size_cowherd，群中的牛数。概念，受孕率。

在此先感谢。

我的代码如下：

size_cowherd<-100 
concep<-0.95 

cows <- function(t=119, mort=0.0005){ 
    G1<- numeric(length = t + 1) 
    G2<- numeric(length = t + 1) 
    G3<- numeric(length = t + 1) 
    G4<- numeric(length = t + 1) 
    G5<- numeric(length = t + 1) 
    G6<- numeric(length = t + 1) 
    G7<- numeric(length = t + 1) 
    G8<- numeric(length = t + 1) 
    G9<- numeric(length = t + 1) 
    MC1<- numeric(length = t + 1) 
    MC2<- numeric(length = t + 1) 
    MC3<- numeric(length = t + 1) 
    MC4<- numeric(length = t + 1) 
    MC5<- numeric(length = t + 1) 
    MC6<- numeric(length = t + 1) 
    MC7<- numeric(length = t + 1) 
    Rest1<- numeric(length = t + 1) 
    Rest2<- numeric(length = t + 1) 
    Rest3<- numeric(length = t + 1) 
    Rest4<- numeric(length = t + 1) 
    Rest5<- numeric(length = t + 1) 
    Rest6<- numeric(length = t + 1) 
    DeadCows <- numeric(length = t + 1) 
    NPreg <- numeric(length = t + 1) 

    G1[1]<- rbinom(1,size_cowherd,concep) 
    G2[1]<- 0 
    G3[1]<- 0 
    G4[1]<- 0 
    G5[1]<- 0 
    G6[1]<- 0 
    G7[1]<- 0 
    G8[1]<- 0 
    G9[1]<- 0 
    MC1[1]<- 0 
    MC2[1]<- 0 
    MC3[1]<- 0 
    MC4[1]<- 0 
    MC5[1]<- 0 
    MC6[1]<- 0 
    MC7[1]<- 0 
    Rest1[1]<-0 
    Rest2[1]<-0 
    Rest3[1]<-0 
    Rest4[1]<-0 
    Rest5[1]<-0 
    Rest6[1]<-0 
    DeadCows[1] <- 0 
    NPreg[1] <- size_cowherd - G1[1] 

    for(step in 1:t){ 
    G2[step+1] <- rbinom(1, G1[step], (1-mort)) 
    G3[step+1] <- rbinom(1, G2[step], (1-mort)) 
    G4[step+1] <- rbinom(1, G3[step], (1-mort)) 
    G5[step+1] <- rbinom(1, G4[step], (1-mort)) 
    G6[step+1] <- rbinom(1, G5[step], (1-mort)) 
    G7[step+1] <- rbinom(1, G6[step], (1-mort)) 
    G8[step+1] <- rbinom(1, G7[step], (1-mort)) 
    G9[step+1] <- rbinom(1, G8[step], (1-mort)) 
    MC1[step+1] <- rbinom(1, G9[step], (1-mort)) 
    MC2[step+1] <- rbinom(1, MC1[step], (1-mort)) 
    MC3[step+1] <- rbinom(1, MC2[step], (1-mort)) 
    MC4[step+1] <- rbinom(1, MC3[step], (1-mort)) 
    MC5[step+1] <- rbinom(1, MC4[step], (1-mort)) 
    MC6[step+1] <- rbinom(1, MC5[step], (1-mort)) 
    MC7[step+1] <- rbinom(1, MC6[step], (1-mort)) 
    Rest1[step+1] <- rbinom(1,MC7[step],(1-mort)) 
    Rest2[step+1] <- rbinom(1,Rest1[step],(1-mort)) 
    Rest3[step+1] <- rbinom(1,Rest2[step],(1-mort)) 
    Rest4[step+1] <- rbinom(1,Rest3[step],(1-mort)) 
    Rest5[step+1] <- rbinom(1,Rest4[step],(1-mort)) 
    Rest6[step+1] <- rbinom(1,Rest5[step],(1-mort)) 
    G1[step+1] <- rbinom(1, Rest6[step], (1-mort)) 
    DeadCows[step+1] <-sum(G1[step]-G2[step+1],G2[step]-G3[step+1],G3[step]- 
    G4[step+1],G4[step]-G5[step+1],G5[step]-G6[step+1],G6[step]- 
    G7[step+1],G7[step]-G8[step+1],G8[step]-G9[step+1],G9[step]- 
    MC1[step+1],MC1[step]-MC2[step+1],MC2[step]-MC3[step+1],MC3[step]- 
    MC4[step+1],MC4[step]-MC5[step+1],MC5[step]-MC6[step+1],MC6[step]- 
    MC7[step+1],MC7[step]-Rest1[step+1],Rest1[step]- 
    Rest2[step+1],Rest2[step]-Rest3[step+1],Rest3[step]- 
    Rest4[step+1],Rest4[step]-Rest5[step+1],Rest5[step]- 
    Rest6[step+1],Rest6[step]-G1[step+1]) 

     if(G1[step]<size_cowherd){ 
     G1[step+1]<- rbinom(1,Rest6[step], concep) 
     NPreg[step+1]<-Rest6[step]-G1[step+1] 
     } 
    } 
    out <-cbind(G1,G2,G3,G4,G5,G6,G7,G8,G9,MC1,MC2,MC3,MC4,MC5,MC6,MC7,Rest1,R 
    est2,Rest3,Rest4,Rest5,Rest6,DeadCows,NPreg) 
    return(out) 
}

下面的输出应该是什么样子的样本。在23个月中，循环再次重启。

G1 G2 G3 G4 G5 G6 G7 G8 G9 MC1 MC2 MC3 MC4 MC5 MC6 MC7 Rest1 Rest2 Rest3 
    Rest4 Rest5 Rest6 DeadCows NPreg 
1 96 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  0  0  0  
     0  0  0  0  4 
2 0 96 0 0 0 0 0 0 0 0 0 0 0 0 0 0  0  0  0  
     0  0  0  0  0 
3 0 0 96 0 0 0 0 0 0 0 0 0 0 0 0 0  0  0  0  
     0  0  0  0  0 
4 0 0 0 96 0 0 0 0 0 0 0 0 0 0 0 0  0  0  0  
     0  0  0  0  0 
5 0 0 0 0 96 0 0 0 0 0 0 0 0 0 0 0  0  0  0  
     0  0  0  0  0 
6 0 0 0 0 0 96 0 0 0 0 0 0 0 0 0 0  0  0  0  
     0  0  0  0  0 
7 0 0 0 0 0 0 96 0 0 0 0 0 0 0 0 0  0  0  0  
     0  0  0  0  0 
8 0 0 0 0 0 0 0 96 0 0 0 0 0 0 0 0  0  0  0  
     0  0  0  0  0 
9 0 0 0 0 0 0 0 0 96 0 0 0 0 0 0 0  0  0  0  
     0  0  0  0  0 
10 0 0 0 0 0 0 0 0 0 96 0 0 0 0 0 0  0  0  0  
     0  0  0  0  0 
11 0 0 0 0 0 0 0 0 0 0 96 0 0 0 0 0  0  0  0  
     0  0  0  0  0 
12 0 0 0 0 0 0 0 0 0 0 0 96 0 0 0 0  0  0  0  
     0  0  0  0  0 
13 0 0 0 0 0 0 0 0 0 0 0 0 96 0 0 0  0  0  0  
     0  0  0  0  0 
14 0 0 0 0 0 0 0 0 0 0 0 0 0 96 0 0  0  0  0  
     0  0  0  0  0 
15 0 0 0 0 0 0 0 0 0 0 0 0 0 0 96 0  0  0  0  
     0  0  0  0  0 
16 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 96  0  0  0  
     0  0  0  0  0 
17 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 96  0  0  
     0  0  0  0  0 
18 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  0 96  0  
     0  0  0  0  0 
19 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  0  0 95  
     0  0  0  1  0 
20 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  0  0  0  
     95  0  0  0  0 
21 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  0  0  0  
     0 95  0  0  0 
22 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  0  0  0  
     0  0 94  1  0

来源

2017-09-26 altfi_SU

您能给出一个预期输出结果的样本吗？ – shea

@shea是的，这里是... –

像这样的东西应该适合你。我认为这里的诀窍是利用矩阵来使簿记更简单一些。

size_cowherd <- 100 

concep <- 0.95 

stage_names <- c(paste0("G",seq(9)), paste0("MC",seq(7)), paste0("Rest",seq(6))) 

cows <- function(size_cowherd, concep, t=220, mort=0.0005, names=stage_names) { 
    n_stages <- length(names) 
    stages <- matrix(0, t, n_stages) 
    dead_cows <- n_preg <- rep(NA, t) 
    stages[1,1] <- rbinom(1, size_cowherd, concep) 
    dead_cows[1] <- 0 
    n_preg[1] <- size_cowherd - stages[1,1] 
    for(tt in 2:t) { 
    stages[tt,1] <- rbinom(1, stages[tt-1,n_stages], 1-mort) 
    for(i in 2:n_stages) { 
     stages[tt,i] <- rbinom(1, stages[tt-1,i-1], 1-mort) 
    } 
    dead_cows[tt] <- sum(stages[tt-1,] - stages[tt,c(2:n_stages,1)]) 
    if(stages[tt-1,1] < size_cowherd) { 
     stages[tt, 1] <- rbinom(1, stages[tt-1,n_stages], concep) 
     n_preg[tt] <- stages[tt-1,n_stages] - stages[tt,1] 
    } 
    } 
    res <- cbind(stages, dead_cows, n_preg) 
    colnames(res) <- c(names, "Dead", "N_Preg") 
    return(res) 
} 

head(cows(100, 0.95), 24) 

     G1 G2 G3 G4 G5 G6 G7 G8 G9 MC1 MC2 MC3 MC4 MC5 MC6 MC7 Rest1 Rest2 Rest3 Rest4 Rest5 Rest6 Dead N_Preg 
    [1,] 96 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0 0  4 
    [2,] 0 96 0 0 0 0 0 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0 0  0 
    [3,] 0 0 96 0 0 0 0 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0 0  0 
    [4,] 0 0 0 96 0 0 0 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0 0  0 
    [5,] 0 0 0 0 96 0 0 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0 0  0 
    [6,] 0 0 0 0 0 96 0 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0 0  0 
    [7,] 0 0 0 0 0 0 95 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0 1  0 
    [8,] 0 0 0 0 0 0 0 95 0 0 0 0 0 0 0 0  0  0  0  0  0  0 0  0 
    [9,] 0 0 0 0 0 0 0 0 95 0 0 0 0 0 0 0  0  0  0  0  0  0 0  0 
[10,] 0 0 0 0 0 0 0 0 0 95 0 0 0 0 0 0  0  0  0  0  0  0 0  0 
[11,] 0 0 0 0 0 0 0 0 0 0 95 0 0 0 0 0  0  0  0  0  0  0 0  0 
[12,] 0 0 0 0 0 0 0 0 0 0 0 95 0 0 0 0  0  0  0  0  0  0 0  0 
[13,] 0 0 0 0 0 0 0 0 0 0 0 0 95 0 0 0  0  0  0  0  0  0 0  0 
[14,] 0 0 0 0 0 0 0 0 0 0 0 0 0 95 0 0  0  0  0  0  0  0 0  0 
[15,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 95 0  0  0  0  0  0  0 0  0 
[16,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 95  0  0  0  0  0  0 0  0 
[17,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 94  0  0  0  0  0 1  0 
[18,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  0 94  0  0  0  0 0  0 
[19,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  0  0 94  0  0  0 0  0 
[20,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  0  0  0 94  0  0 0  0 
[21,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  0  0  0  0 94  0 0  0 
[22,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  0  0  0  0  0 94 0  0 
[23,] 92 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0 0  2 
[24,] 0 92 0 0 0 0 0 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0 0  0

来源

2017-09-28 23:09:35

谢谢你，马克。这非常有帮助。 @mark s –

如何使用这种方法来生成更简单和可读的输出，以便列名称代表阶段序列的总和，例如， PregCows会在相同的列名称下显示九行的模拟等等？ –

@Francis_SU你能否澄清你的要求？例如，你是否想要一个有10行和22列的表格，其中每行是“运行”，还是你想要一个有1行和22列的表格，其中每列是所有列的总和？ –

R中的长时间序列函数的简化/浓缩

回答

相关问题