[R-SIG-Finance] objective function optimization in JAGS/BUGS rather then in optim-like R tools

[Alex Bird]

Hi there!

  I'm not sure whether this place is the right one to ask my question
so sorry in advance.
  My problem consists of two steps: on the first one I fit some model
parameters using R2jags and do n-days ahead prediction, then on the
second step I use doptim package to optimize my objective function
based on the predicted data. So I wonder is there a way to find
minimum/maximum of some function directly in JAGS or Win/OpenBUGS.

  To get an idea of what I am looking for please take a look at the
following code

require(R2jags)

### first step - fitting params
#########################################
#helper function
write.jags.model<-function(model,filename="model.bug",dir=getwd()){
    old.dir <- getwd()
    setwd(dir)
    on.exit(setwd(old.dir))
    if (file.exists(filename)) {
        sn <- unlist(strsplit(filename, "\\."))
        if (length(sn) > 2) {
            sn[(length(sn) - 1)] <- paste(sn[-length(sn)], collapse = ".")
            sn <- sn[-c(1:(length(sn) - 2))]
        }
        ff <- tempfile("model", "")
        filename2 <- paste(substr(ff, 2, nchar(ff)), "bug", sep = ".")
        if (file.exists(filename2)) {
            while (!file.exists(filename2)) {
                ff <- tempfile("model", "")
                filename2 <- paste(substr(ff, 2, nchar(ff)),
                  "bug", sep = ".")
            }
        }
    }
    else {
        filename2 <- filename
    }
    if (inherits(model, "custommodel")) {
        writeLines(model, filename2)
    }
    else {
        write.model(model, filename2)
    }
    invisible(filename2)
}

#simulated data
N<-100
y<-rnorm(N,.001,.01)

#model to fit (actually wrong one)
jmodel<-function(){
	for (i in 1:N){
		size[i]~dnorm(muj,precj)
		j[i]~dbern(lambda)
		dummy[i]<-mu+j[i]*size[i]
		y[i]~dnorm(dummy[i],prec)
	}
	mu~dnorm(0,0.01)
	prec~dgamma(2,0.01)
	lambda~dbeta(1,5)
	muj~dnorm(0,0.01)
	precj~dgamma(5,20)
	
	for (i in 1:5){
		size.pred[i]~dnorm(muj,precj)
		j.pred[i]~dbern(lambda)
		dummy.pred[i]<-mu+j.pred[i]*size.pred[i]
		y.pred[i]~dnorm(dummy.pred[i],prec)		
	}
}
model.file<-write.jags.model(jmodel)
jparams<-c('mu','prec','lambda','muj','precj','y.pred')
jags.data<-list('y','N')
jfit<-jags(data=jags.data,inits=NULL,jparams,n.iter=10000,model.file=model.file,DIC=T,n.chains=3)


### second step - function optimization based on y.pred
#########################################
y.pred<-colMeans(jfit$BUGSoutput$sims.list$y.pred)
x.pred<-rnorm(5,.001,.01)

#some objective function - just to get an idea
utf<-function(w=1/2){-c(mean(y.pred),mean(x.known))%*%c(w,1-w)}
w.opt<-optimize(utf,lower=-1,upper=1)$minimum
w.opt


where the second step I want to put somehow directly into jmodel

Thanks!
Alex