我希望建立一個人口動態模型，其中每個引數值對應于當天的溫度。例如
簡單模型
library(deSolve)
set.seed(1)
pars <- c(alpha = 1, beta = 0.2, delta = 0.5, gamma = 0.2)
lv_model <- function(pars, times = seq(0, 50, by = 1)) {
# initial state
state <- c(x = 1, y = 2)
# derivative
deriv <- function(t, state, pars) {
with(as.list(c(state, pars)), {
d_x <- alpha * x - beta * x * y
d_y <- delta * beta * x * y - gamma * y
return(list(c(x = d_x, y = d_y)))
})
}
# solve
ode(y = state, times = times, func = deriv, parms = pars)
}
lv_results <- lv_model(pars = pars, times = seq(0, 50, by = 1))
我現在想使用一系列每日溫度
DailyTemperature<-floor(runif(50,0,40))
并使引數值成為溫度的函式
TraitTemperature<-seq(1,40,1)
#trait responses to temperature
alpha<- abs(rnorm(40,mean = 0.5,sd=1))
beta<- abs(rnorm(40,mean = 0.2,sd=0.5))
delta<-abs(rnorm(40,mean=1,sd=2))
gamma<- seq(0.025,1,0.025)
parameters<-as.data.frame(cbind(TraitTemperature,alpha,beta,delta,gamma))
這樣對于迭代的每個時間步，它都會查看每日溫度，然后在引數資料框中找到相應的溫度值。
回顧檔案，我看到了if/else在特定時間步和強制函式的使用中想要更改單個引數時使用的陳述句，但我認為它們不適用于這里。
我希望這是有道理的，我對如何使它起作用的想法很感興趣。到目前為止，我還嘗試使用 afor loop來遍歷每日溫度串列，然后使用match函式來識別值，但這并未涉及每日時間步長。
uj5u.com熱心網友回復：
這是一種可能的方法，使用DailyTemperature強制，然后使用parameters資料框作為查找表。match這里不需要，只要溫度是整數，并且資料幀中的溫度對應資料幀的行索引即可。
我想補充一點，原則上，不連續強迫會使模型變慢，因為 ODE 根據定義是連續的。幸運的是，由于求解器非常健壯，它應該適用于實際應用：
library(deSolve)
set.seed(1)
deriv <- function(t, state, pars) {
pars <- parameters[DailyTemperature[floor(t 1)], 2:5]
#print(pars)
with(as.list(c(state, pars)), {
d_x <- alpha * x - beta * x * y
d_y <- delta * beta * x * y - gamma * y
list(c(x = d_x, y = d_y), alpha=alpha, beta=beta, gamma=gamma, delta=delta)
})
}
state <- c(x = 1, y = 2)
times = seq(0, 50, by = 1)
# pars <- c(alpha = 1, beta = 0.2, delta = 0.5, gamma = 0.2)
parameters <- data.frame(
TraitTemperature = seq(1,40,1),
alpha = abs(rnorm(40,mean = 0.5,sd=1)),
beta = abs(rnorm(40,mean = 0.2,sd=0.5)),
delta = abs(rnorm(40,mean=1,sd=2)),
gamma = seq(0.025,1,0.025)
)
DailyTemperature <- floor(runif(51, 0, 40)) # one more because start zero
out <- ode(y = state, times = times, func = deriv, parms = pars)
plot(out)
引數作為串列
在上面的例子中，pars傳入的變數ode只是被pars從全域變數派生的parameters和覆寫了DailyTemperature。這是可行的，但也可以考慮將兩者作為串列傳遞給deriv函式。
deriv <- function(t, state, p) {
parameters <- p[[1]]
DailyTemperature <- p[[2]]
parms <- parameters[DailyTemperature[floor(t 1)], 2:5]
# ...
}
接著：
out <- ode(y = state, times = times, func = deriv,
parms = list(parameters, DailyTemperature))
uj5u.com熱心網友回復：
這似乎適用于使用當前可重現代碼的索引：
set.seed(1)
deriv <- function(t, state, pars) {
pars<- parameters[match(parameters$TraitTemperature[parameters[2:5]],DailyTemperature),]
#print(pars)
with(as.list(c(state, pars)), {
d_x <- alpha * x - beta * x * y
d_y <- delta * beta * x * y - gamma * y
list(c(x = d_x, y = d_y), alpha=alpha, beta=beta, gamma=gamma, delta=delta)
})
}
state <- c(x = 1, y = 2)
times = seq(0, 50, by = 1)
# pars <- c(alpha = 1, beta = 0.2, delta = 0.5, gamma = 0.2)
parameters <- data.frame(
TraitTemperature = seq(1,40,1),
alpha = abs(rnorm(40,mean = 0.5,sd=1)),
beta = abs(rnorm(40,mean = 0.2,sd=0.5)),
delta = abs(rnorm(40,mean=1,sd=2)),
gamma = seq(0.025,1,0.025)
)
DailyTemperature <- floor(runif(51, 0, 40)) # one more because start zero
out <- ode(y = state, times = times, func = deriv, parms = pars)
plot(out)
但是如果我增加值的解析度，例如
# Parameter datasets
parameters <- data.frame(
TraitTemperature = seq(0.1,40,0.1),
alpha = abs(rnorm(400,mean = 0.5,sd=1)),
beta = abs(rnorm(400,mean = 0.2,sd=0.5)),
delta = abs(rnorm(400,mean=1,sd=2)),
gamma = seq(0.0025,1,0.0025)
)
# random daily temperature dataset
DailyTemperature <- round(runif(51, 0, 40),1) # one more because start zero
然后我在某些時間步驟得到 NA
uj5u.com熱心網友回復：
OP 擴展了她/他的問題，因此這可能是開始新執行緒的首選方式，但為了提供快速反饋，我嘗試給出另一個答案。
有幾種方法可用于在 2D（時間和溫度）中進行插值或索引。我的首選方法是創建一個 2D 模型，然后使用 2D 插值方法。如果引數表面是光滑的，而不僅僅是給定示例中的隨機，則此方法效果最佳。但是，為了簡單起見，還可以使用舍入和查表。如果值不是整數，由于舍入效應和有限的精度（IEEE 數字格式），精確比較通常不起作用，因此可以使用which.min：
DailyTemperature <- round(runif(51, 0, 40), 1)
TraitTemperature <- seq(0, 40, by=0.1)
N <- length(TraitTemperature)
parameters <- data.frame(
TraitTemperature = TraitTemperature,
alpha = abs(rnorm(N, mean = 0.5, sd=1)),
beta = abs(rnorm(N, mean = 0.2, sd=0.5)),
delta = abs(rnorm(N, mean=1,sd=2)),
gamma = seq(0.025, 1, length.out=N)
)
t <- 17
actualTemp <- DailyTemperature[floor(t 1)]
actualTemp
pars <- parameters[which.min(abs(actualTemp - parameters$TraitTemperature)), 1:5]
head(pars)
轉載請註明出處，本文鏈接：https://www.uj5u.com/qukuanlian/388857.html
上一篇：如何使它成為嵌套的For回圈？