R Параллельный EasyAB C или smfsb: AB C -SMC - PullRequest
1 голос
/ 19 июня 2020

У меня есть простой ODE с некоторыми неизвестными параметрами (r, C, d и g), который предсказывает разрушение и повторный рост бактерий на поверхности после очистки. У меня есть экспериментальные данные о количестве бактерий в различные моменты времени (0, 1, 2, 4, 8 и 24 часа) после очистки. Я пытаюсь использовать пакет EasyAB C для оценки распределения параметров r, C, d и g с помощью подхода приближенного байесовского вычисления последовательного моделирования Монте-Карло.

enter image description here

Установить необходимые пакеты

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 и во многих отношениях это более простой пример, чем приведенный в документации, но я не могу понять, как передать в него свою функцию. У кого-нибудь есть подобный пример?

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...