Dear All,
Currently I am running the following code:
library(stats4)
library(odesolve)
library(rgenoud)
Input<-data.frame(SUB=c(1),time=c(0.5,3,10,15),lev=c(2.05,12.08,9.02,8))
XD<-500
IT<-3
diffeqfun<-function(time, y, parms) {
if(time<=IT)
dCpdt <- (XD/IT)/parms["Vol"] -
(parms["Vm"]/parms["Vol"])*y[1]/(parms["Km"]/parms["Vol"]+y[1])
else
dCpdt <-
-(parms["Vm"]/parms["Vol"])*y[1]/(parms["Km"]/parms["Vol"]+y[1])
list(dCpdt)
}
modelfunc<-function(time,Vm,Km,Vol) {
out <-
lsoda(0,c(0,time),diffeqfun,parms=c(Vm=Vm,Km=Km,Vol=Vol),rtol=1e-5,atol=1e-5)
out[-1,2]
}
objectfunc <- function(par) {
out <- modelfunc(Input$time, par[1], par[2], par[3])
gift <- which( Input$lev != 0 )
sum((Input$lev[gift]-out[gift])^2)
}
gen <-
genoud(objectfunc,nvars=3,max=FALSE,pop.size=10,max.generations=100,wait.generations=100,
starting.value=c(40,8,12),BFGS=FALSE,
print.level=0,boundary.enforcement=0,Domains=matrix(c(0,0,0,100,100,100),3,2),MemoryMatrix=TRUE)
outgen<-c(gen$par[1],gen$par[2],gen$par[3])
opt<-optim(c(gen$par[1],gen$par[2],gen$par[3]),objectfunc,method="Nelder-Mead")
outopt<-c(opt$par[1],opt$par[2],opt$par[3])
fm<-nls(lev~modelfunc(time,Vm,Km,Vol),data=Input,start=list(Vm=opt$par[1],Km=opt$par[2],Vol=opt$par[3]),trace=TRUE,nls.control(tol=1))
x<-Input$time
y<-Input$lev
cal<-predict(fm,list(time=x))
plot(x,y,type="l",col="green")
lines(x,cal,col="red")
as you can see, the parameter optimization using this code is done via a
non-linear approach using the Nelder-Mead optimization. I was wondering if you
could give me some ideas on how to rewrite this code so that it could be run
using OpenBUGS via R2OpenBUGS? The issue I am hung up on is how to write
the differential equation into a code that could run through R2OpenBUGS and
OpenBUGS. Below is an example for an analytical solution for a different model I
currently use that is run through R2OpenBUGS and OpenBUGS that generates
excelent results (to be able to run this code OpenBUGS installed, also please
note this is only an example, I really need to reproduce the ideas from above
using the bayesian approach):
library(R2OpenBUGS)
D <-400
tin <-1
VDPRE <-0.5
LWEIGHT <-68
CRCL <-60
tau <-12
T <-2
ts <-c(2,8)
c <-c(15,6)
mfsssc <- function() {
for (i in 1:1) {
for (j in 1:T) {
c[j]~dnorm(mu[j],100)
mu[j]<-((D/tin)*(1/cl))*(((1-exp(-(cl/v)*tin))*exp(-(cl/v)*(ts[j])))/(1-exp(-(cl/v)*tau)))
}
cla~dnorm(0, 6.25)
va~dnorm(0, 6.25)
theta1<-0.615*(CRCL*0.06)+0.7194
theta2<-VDPRE*LWEIGHT
cl<-exp(cla)*theta1
v<-exp(va)*theta2
}
}
res <- bugs(list(T = T, c = c, tau = tau, ts = ts, tin = tin, D = D, CRCL =
CRCL, VDPRE = VDPRE, LWEIGHT = LWEIGHT), model.file = mfsssc, inits = NULL,
parameters.to.save=c("v", "cl"), n.chains=2, n.iter=10000)
finalsssc <- res$summary
I hope my question made enough sense, but if not please ask and I will try to
further explain my goals here. All your input would be greatly apreciated,
Sincerely,
Andras
[[alternative HTML version deleted]]