Messagepar Jue Binette » 10 Sep 2008, 14:30
J'essais tout ça... en attendant, voici le script auquel la fonction wbs refère
wbs<-function(data0=thewho, model = '1diff2stateM.bug', samples = 2000, burn = 1000, thin = 2, sim = FALSE, sim.data = NULL, dim = 2, debug = FALSE, codaPkg=F, chains = 2, switch = 2){
### List of proceedures to prepare data and initial values for a number of parameters.
loc.tmp<-data0
oidx<-max(loc.tmp$Oidx)
Ridx<-loc.tmp$Ridx
RegN<-max(Ridx)
RegO<-RegN+1
numtracks <- length(loc.tmp$Ridx)-1
y<-loc.tmp$y
row <- unique(1:nrow(y))
na.row <- row[is.na(y)]
na.row <- na.row[1:(length(na.row)/2)]
# y.init <- y
# y.init[!is.na(y.init)] <- 9999
# tmp <- sapply(1:length(na.row), function(i){c(mean(y[c(na.row[i] - 1, na.row[i] + 1),1], na.rm = T), mean(y[c(na.row[i] - 1, na.row[i] + 1),2], na.rm = T))})
# tmp[tmp=='NaN'] <- rep(c(mean(tmp[1,], na.rm = T), mean(tmp[2,], na.rm = T)), length(tmp[tmp=='NaN'])/2)
# y.init[na.row,] <- t(tmp)
# y.init[y.init == 9999] <- NA
bmode.init <- unlist(lapply(1:numtracks, function(i){c(rep(1, (Ridx[i+1] - Ridx[i] - 1)), NA)}))
bmode.init <- bmode.init[-length(bmode.init)]
tmp<-strsplit(deparse(substitute(data0)),"\\.")[[1]][1]
### Final preparation of initials, data, and paramaters to follow through the MCMC proceedure
data1 <- list(y = loc.tmp$y, idx = loc.tmp$idx, j = loc.tmp$j, Ridx = loc.tmp$Ridx, numtracks = numtracks, itau = loc.tmp$itau, nu = loc.tmp$nu, Omega = matrix(c(1,0,0,1),2,2))
inits <- list(list(iSigma = matrix(c(1, 0, 0, 1), 2, 2), alpha = c(0.9,0.1), vfact= rep(0.5,numtracks), lambda = matrix(rep(c(0.5, NA),numtracks),numtracks,2,byrow=T), gamma1 = c(0.5,0.5), tprior = c(0.5,0.5), bmode = bmode.init, x = matrix(rep(loc.tmp$y[1,],RegO-2),RegO-2,2,byrow=T)), list(iSigma = matrix(c(1, 0, 0, 1), 2, 2), alpha = c(0.9,0.1), vfact= rep(0.5,numtracks), lambda = matrix(rep(c(0.5, NA),numtracks),numtracks,2,byrow=T), gamma1 = c(0.5,0.5), tprior = c(0.5,0.5), bmode = bmode.init, x = matrix(rep(loc.tmp$y[1,],RegO-2),RegO-2,2,byrow=T))) #bmode = bmode.init,
parameters <- c('Sigma', 'tprior', 'gamma', 'bmode', 'alpha', 'x', 'vfact')
### End final preparation of initials, data, and parameters to follow.
ptm<-proc.time() # stopwatch of CPU time
### Call WinBUGS
test1<-bugs(data1, inits, parameters.to.save=parameters, model, n.chains = chains, n.iter = samples, n.burnin = burn, n.thin = thin, debug = debug, digits = 8, codaPkg = codaPkg, bugs.directory = "C:/Documents and Settings/harveyv/Mes documents/WinBUGS14/")
wbtime<-proc.time()-ptm # stopwatch of CPU time
### Move results to R workspace
test1<-list(summary = test1$summary, median = test1$median, mean=test1$mean, DIC = test1$DIC, pD = test1$pD, sims.matrix = as.data.frame(test1$sims.matrix), wbtime = wbtime[3], tstep=loc.tmp$timestep, btfac=loc.tmp$btfac, bt=loc.tmp$bt, model=model, Ridx = data0$Ridx, rmonth = data0$rmonth, lc=data0$lc, sead = data0$sead, seallist = data0$seallist, Oidx = data0$Oidx, POidx = data0$POidx)
assign(paste(tmp,'.wb',sep=''),test1,pos=1)
print(test1$wbtime,digits.summary=3)
}
encore une fois merci!