У меня есть простой ODE с некоторыми неизвестными параметрами (r, C, d и g), который предсказывает разрушение и повторный рост бактерий на поверхности после очистки. У меня есть экспериментальные данные о количестве бактерий в различные моменты времени (0, 1, 2, 4, 8 и 24 часа) после очистки. Я пытаюсь использовать пакет EasyAB C для оценки распределения параметров r, C, d и g с помощью подхода приближенного байесовского вычисления последовательного моделирования Монте-Карло.
Установить необходимые пакеты
if("deSolve" %in% rownames(installed.packages())==FALSE){install.packages("deSolve"); require(deSolve)}else{require(deSolve)}
Функция уравнения:
ode_model <- function(time_space, initial_contamination, parameters){
with(
as.list(c(initial_contamination, parameters)),{
dContamination <- Contamination*r*(1-Contamination/C)-d*exp(-g*time_space)*Contamination
return(list(dContamination))
}
)
}
Подсчет бактерий в 5 временных точках
deterministic_run<-function(precision,initial_contamination,parameters){
tmax = 24
time_space = seq(0, tmax, length.out = precision+1)
sim<- ode(initial_contamination, time_space, ode_model, parameters, method = "radau",atol = 1e-4, rtol = 1e-4)
num_at_1=sim[(precision*1/50.0),2]
num_at_2=sim[(precision*2/50.0),2]
num_at_4=sim[(precision*4/50.0),2]
num_at_8=sim[(precision*8/50.0),2]
num_at_24=sim[(precision*24/50.0),2]
return(cbind(num_at_1,num_at_2,num_at_4,num_at_8,num_at_24))
}
Это запускает решатель ODE
trial_r = runif(1,1E-1,4)
trial_C = runif(1,6,15)
trial_d = runif(1,1,4)
trial_g = runif(1,1E-2,1)
trial_l = runif(1,1E-5,2)
parameters=c(r=trial_r,C=trial_C,d=trial_d,g=trial_g,l=trial_l)
initial_contamination=59
precision=5000
one_run = deterministic_run(precision,initial_contamination,parameters)#
Plot
plot(one_run)
Я пытаюсь использовать пакет EasyAB C: http://easyabc.r-forge.r-project.org для оценки параметров в ODE и во многих отношениях это более простой пример, чем приведенный в документации, но я не могу понять, как передать в него свою функцию. У кого-нибудь есть подобный пример?