deSolve包的參數可以包含一個矩陣嗎？

我正在嘗試編碼年齡分層的SEIR模型;也就是說，在我的微分方程中，我有一個質量作用的參數，它是在20個年齡段中beta *（受感染比例）*（易感人數）的總和。透射係數（β）由接觸矩陣計算。聯繫矩陣有20列和行，它們代表年齡等級（行=人i，列=人j），幷包含任何年齡段中兩個人之間接觸的概率。我設計了它並將其讀入R.我的問題是我不知道如何（或如果）我可以在deSolve中使用我的參數內的矩陣。此代碼下面我寫的不工作，我認爲，由於矩陣的/我得到這個錯誤：deSolve包的參數可以包含一個矩陣嗎？

Error in beta * S : non-numeric argument to binary operator

之前，我與它傻瓜太多了，我想知道是否可以使用這樣的矩陣作爲這個模型的一個參數。

mat <-as.matrix(read.csv("H:/IBS 796R/contactmatrix.csv", header=F)) 

times <- seq(0,20,by=1/52) 
parameters <- c(mu=0,v=1/75,N=1,p=0,delta=2.4,beta=mat*0.04,sigma=1/8,gamma=1/15) 
xstart <- c(S=0.06,E=0,I=0.001,R=0) 

SEIR0 <- function(t,x,parameters){ 
    S=x[1] 
    E=x[2] 
    I=x[3] 
    R=x[4] 
    with(as.list(parameters), { 
     dS=v*S -beta*S*I/N -delta*S 
     dE=beta*S*1/N -E*(sigma+delta) 
     dI=sigma*E -I*(gamma+delta) 
     dR=gamma*I-delta*R 
     res=c(dS,dE,dI,dR) 
     list(res) 
    }) 
} 

out <- as.data.frame(lsoda(xstart,times,SEIR0,parameters))

另外，如果我打印參數，這是什麼看起來像公測：

$beta.V1 
[1] 4e-04 4e-04 4e-04 4e-04 4e-04 4e-04 4e-04 4e-04 4e-04 4e-04 4e-04 
[12] 4e-04 4e-04 8e-03 8e-03 8e-03 8e-03 8e-03 8e-03 8e-03 

$beta.V2 
[1] 4e-04 4e-04 4e-04 4e-04 4e-04 4e-04 4e-04 4e-04 4e-04 4e-04 4e-04 
[12] 4e-04 4e-04 8e-03 8e-03 8e-03 8e-03 8e-03 8e-03 8e-03

....到$ beta.V20。所以我認爲它創建了20個向量，每個向量有20個參數...我認爲每個原始矩陣「mat」乘以常數0.04是一排？但是，當我將mat * 0.04乘以「參數」外，我得到預期的矩陣。我正在爲如何使用deSolve實現這些方程而苦苦掙扎，並希望得到關於是否可能的任何建議。提前致謝。

來源

2013-03-11 user2100907

錯誤發生在這條線：

dS=v*S -beta*S*I/N -delta*S

non-numeric argument to binary operator意味着你試圖通過一個數字乘以例如功能。您可以通過I*1

Error in I * 1 : non-numeric argument to binary operator`

在這裏重現，R可以」找到測試版，和β被解釋爲數學的特殊功能，所以錯誤。您需要定義參數

# a list 
list(mu=0,v=1/75,N=1,p=0,delta=2.4,beta=mat*0.04,sigma=1/8,gamma=1/15)

和

## you get a vector mu,N,p,delta,beta1,bet2,... 
c(mu=0,v=1/75,N=1,p=0,delta=2.4,beta=mat*0.04,sigma=1/8,gamma=1/15)

我想你甚至可以重寫你的函數爲：

SEIR0 <- function(t,x,parameters){ 
    with(as.list(c(parameters, x)), { 
    dS = v*S -beta*S*I/N -delta*S ## matrix 
    dE = beta*S*1/N -E*(sigma+delta) ## matrix 
    dI = sigma*E -I*(gamma+delta) 
    dR = gamma*I-delta*R 
    res = c(dS,dE,dI,dR) 
    list(res)      ## different of the structure of xstart 
    }) 
}

這將解決上述問題，但ODE不會因爲SEIR0返回的衍生物的數量必須等於初始條件的長度xstart矢量（這裏是4）。

我建議，例如：

res <- c(dS=mean(dS),dE=mean(dE),dI=dI,dR=dR) 
    list(res)

來源

2013-03-11 06:01:17 agstudy

deSolve包的參數可以包含一個矩陣嗎？

回答

相關問題