Menu
▾
▴
[Bugs-r-devel] SF.net SVN: bugs-r: [3] trunk/R2WinBUGS
|
From: <gg...@us...> - 2006-11-24 13:36:30
|
Revision: 3
http://svn.sourceforge.net/bugs-r/?rev=3&view=rev
Author: ggorjan
Date: 2006-11-24 05:24:51 -0800 (Fri, 24 Nov 2006)
Log Message:
-----------
populating SVN for R2WinBUGS
Added Paths:
-----------
trunk/R2WinBUGS/Changes
trunk/R2WinBUGS/Changes0
trunk/R2WinBUGS/DESCRIPTION
trunk/R2WinBUGS/NAMESPACE
trunk/R2WinBUGS/R/
trunk/R2WinBUGS/R/as.bugs.array.R
trunk/R2WinBUGS/R/attach.all.R
trunk/R2WinBUGS/R/bugs.R
trunk/R2WinBUGS/R/bugs.data.R
trunk/R2WinBUGS/R/bugs.inits.R
trunk/R2WinBUGS/R/bugs.log.R
trunk/R2WinBUGS/R/bugs.plot.inferences.R
trunk/R2WinBUGS/R/bugs.plot.summary.R
trunk/R2WinBUGS/R/bugs.run.R
trunk/R2WinBUGS/R/bugs.script.R
trunk/R2WinBUGS/R/bugs.sims.R
trunk/R2WinBUGS/R/bugs.update.settings.R
trunk/R2WinBUGS/R/conv.par.R
trunk/R2WinBUGS/R/decode.parameter.name.R
trunk/R2WinBUGS/R/formatdata.R
trunk/R2WinBUGS/R/monitor.R
trunk/R2WinBUGS/R/openbugs.R
trunk/R2WinBUGS/R/plot.bugs.R
trunk/R2WinBUGS/R/print.bugs.R
trunk/R2WinBUGS/R/read.bugs.R
trunk/R2WinBUGS/R/wineutils.R
trunk/R2WinBUGS/R/write.datafile.R
trunk/R2WinBUGS/R/write.model.R
trunk/R2WinBUGS/data/
trunk/R2WinBUGS/data/schools.txt
trunk/R2WinBUGS/inst/
trunk/R2WinBUGS/inst/CITATION
trunk/R2WinBUGS/inst/doc/
trunk/R2WinBUGS/inst/doc/R2WinBUGS.Rnw
trunk/R2WinBUGS/inst/doc/R2WinBUGS.pdf
trunk/R2WinBUGS/inst/doc/RdRW.sty
trunk/R2WinBUGS/inst/doc/Z.cls
trunk/R2WinBUGS/inst/doc/benzolsw.pdf
trunk/R2WinBUGS/inst/doc/bugs.tex
trunk/R2WinBUGS/inst/doc/countssw.pdf
trunk/R2WinBUGS/inst/doc/expectedsw.pdf
trunk/R2WinBUGS/inst/doc/literatur.bib
trunk/R2WinBUGS/inst/doc/literatur.bst
trunk/R2WinBUGS/inst/doc/plot.pdf
trunk/R2WinBUGS/inst/model/
trunk/R2WinBUGS/inst/model/schools.txt
trunk/R2WinBUGS/man/
trunk/R2WinBUGS/man/as.bugs.array.Rd
trunk/R2WinBUGS/man/attach.all.Rd
trunk/R2WinBUGS/man/bugs.Rd
trunk/R2WinBUGS/man/bugs.data.Rd
trunk/R2WinBUGS/man/bugs.inits.Rd
trunk/R2WinBUGS/man/bugs.log.Rd
trunk/R2WinBUGS/man/bugs.plot.Rd
trunk/R2WinBUGS/man/bugs.run.Rd
trunk/R2WinBUGS/man/bugs.script.Rd
trunk/R2WinBUGS/man/bugs.sims.Rd
trunk/R2WinBUGS/man/bugs.update.settings.Rd
trunk/R2WinBUGS/man/decode.parameter.name.Rd
trunk/R2WinBUGS/man/monitor.Rd
trunk/R2WinBUGS/man/openbugs.Rd
trunk/R2WinBUGS/man/plot.bugs.Rd
trunk/R2WinBUGS/man/print.bugs.Rd
trunk/R2WinBUGS/man/read.bugs.Rd
trunk/R2WinBUGS/man/schools.Rd
trunk/R2WinBUGS/man/write.datafile.Rd
trunk/R2WinBUGS/man/write.model.Rd
Added: trunk/R2WinBUGS/Changes
===================================================================
--- trunk/R2WinBUGS/Changes (rev 0)
+++ trunk/R2WinBUGS/Changes 2006-11-24 13:24:51 UTC (rev 3)
@@ -0,0 +1,71 @@
+Changes to R2WinBUGS:
+=====================
+
+During the process of packaging R2WinBUGS_0.1 - R2WinBUGS_0.2-1,
+quite a lot of changes had been made. Those changes are not
+documented anywhere ...
+
+
+Update 0.2-2 (28 Apr 2004):
+- schools data: original (see references) instead of the rounded data
+
+Update 0.2-3 (10 Sept 2004):
+- bugs.data.inits split to bugs.inits and bugs.data, the latter exported
+ from the Namespace.
+ Now we can use already written data files in bugs().
+
+Update 0.2-4 (05 Oct 2004):
+- bugs.script() did not work for large n.iter values in update step
+ (no scientific notation allowed)
+
+Update 0.2-5 (20 Oct 2004):
+- bugs() and bugs.script() have a new argument bin that allows to specify
+ a number of iterations. After each "bin" iterations the coda files are saved.
+
+Update 0.2-6 (18 May 2005):
+- bugs() changes in order to return a file names of coda output files
+- new read.bugs() returns a coda mcmc.list object, if codaPkg=TRUE.
+
+Update 0.2-8 (30 May 2005):
+- bugs passed DIC to bugs.script in order to be able to disable it
+
+Update 0.2-9 (26 July 2005):
+- bugs has new argument clearWD
+
+Update 1.0-0 (05 Aug 2005):
+- attach.all(), detach.all(), attach.bugs() and detach.bugs() added/changed
+ more or less according to Andrew Gelman's current bugs.R
+
+Update 1.0-1 (14 Nov 2005):
+- make inits=NULL work (again ?)
+
+Update 1.1-0 (14 Dec 2005):
+- Contribution by Ben Bolker and Yun Yan's rbugs package:
+ make R2WinBUGS work under WINE
+
+Update 1.1-1 (17 Feb 2006):
+- WINE tweaks (mainly by Gregor Gorjanc)
+
+Update 2.0-0 (08 May 2006):
+- bugs() doc fix/updates on scrambling
+- bugs.run() has new arg useWINE (by Gregor Gorjanc)
+- bugs() and bugs.script() patched for WINEPATH issue
+- bugs.script() changed to save log file in ASCII
+- new function bugs.log() by Gregor Gorjanc
+- new functions as.bugs.array, openbugs and sort.name by
+ Jouni Kerman and Andrew Gelman
+- new function write.model() based on ideas from Jouni Kerman
+
+Update 2.0-1 (26 May 2006):
+- some wine patches for 2.2-0 by Gregor Gorjanc
+
+Update 2.0-2 (26 July 2006):
+- changes for DIC, making use of BUGS internal calculations
+- some doc fixes
+
+Update 2.0-3 (06 October 2006):
+- \\. -> [.] in regular expressions
+
+Update 2.0-4 (01 November 2006):
+- print.bugs / plot.bugs documentation fixes
+- write.model() fix
Added: trunk/R2WinBUGS/Changes0
===================================================================
--- trunk/R2WinBUGS/Changes0 (rev 0)
+++ trunk/R2WinBUGS/Changes0 2006-11-24 13:24:51 UTC (rev 3)
@@ -0,0 +1,142 @@
+During the process of packaging R2WinBUGS_0.1, quite a lot of changes had
+been made. Those changes are not documented anywhere ...
+
+
+Changes prior to R2WinBUGS_0.1:
+===============================
+
+Update 30 Oct 2003:
+ 1. Minor change to a return() statement to be compatible with R 1.8
+ 2. Just a warning (from Ben Goodrich): if you are running Bugs.R inside
+ a loop, you have to be careful with the data() and the inits()
+ functions since they will not necessarily recognize locally-defined
+ variables. One workaround is to define the variables used in data()
+ and inits() using global assignments (<<-), but this can sometimes
+ make the program run slower.
+Update 29 Aug 2003:
+ 1. Fixed "bugs.data.inits" function so you can use data that have the
+ same names as R functions.
+ 2. Changed T and F to TRUE and FALSE everywhere in case the variables
+ T and F are used as data in the main program
+ 3. Caution: if you are entering the data as a list of variable names
+ (see 10 Apr 2003 update, item 1), the data to be input into must
+ be global variables. This can be relevant if you are running bugs()
+ inside an R function.
+ 4. Caution: bugs() has difficulty processing ragged arrays. It is
+ better to save a whole matrix (e.g., "theta") rather than parts
+ (e.g., "theta[1:10,1]", "theta[1:5,2]"). If you want to save
+ part of a vector, you should do it as "theta[1:2]", not "theta[1]",
+ "theta[2]".
+Update 30 Apr 2003: added time monitoring
+Update 29 Apr 2003:
+ 1. The "attach.all" function (no longer called "attach2") overwrites
+ so that all components of a list or data frame get attached.
+ 2. Program now looks in the directory /winbug~1/ rather than /winbug~2/
+ 3. Graphics parameters for margins are returned to their original state
+ at the end of the program.
+ 4. Added "digits.summary" option to the numerical display.
+ 5. Added "last.values" output: a list that can be input as "inits"
+ if you want to run the simulations longer.
+Update 13 Apr 2003: fixed new bug in round.bugs(). Now all numbers are
+ saved in scientific notation.
+Update 10 Apr 2003:
+ 1. It is now possible to enter the data as a list of variable names.
+ For example, before you had to enter data as,
+ data <- list (n=8, y=c(28,8,-3,7,-1,1,18,12))
+ or
+ n <- 8
+ y <- c(28,8,-3,7,-1,1,18,12)
+ data <- list (n=n, y=y)
+ Now you can enter the data as,
+ n <- 8
+ y <- c(28,8,-3,7,-1,1,18,12)
+ data <- list ("n", "y")
+ The bugs() function will figure out which method you are using (based
+ on whether "data" is a list of numbers or a vector of character
+ strings).
+ This doesn't look like much, but it's convenient when you're entering
+ lots of data!
+ 2. It is now possible to enter the initial values as a function,
+ so as to automatically a random list of inits for each of the chains.
+ For example, in the 8-schools example below, we can do:
+ inits <- function()
+ list (theta=rnorm(J,0,1), mu.theta=rnorm(1,0,100), sigma.theta=runif(1,0,100))
+
+ to set up the inits as a function (rather than setting up n.chains
+ sets of specific initial values). Then, the function call,
+ schools.sim <- bugs (data, inits, parameters, "schools.txt", n.chains=3, n.iter=1000)
+ automatically sets up 3 sets of initial values (each a list of
+ theta, mu.theta, sigma.theta).
+ 3. Bug in the initial rounding (the function round.bugs()) has been fixed.
+ Thanks for Mark Clements for finding the bug and fixing it!
+ Also, we have set the default rounding to 5 digits instead of 2.
+Update 01 Apr 2003: use layout() rather than split.screen() for graphical
+ display
+Update 18 Mar 2003:
+ 1. Get the Bugs configuration information from the original file
+ (Registry_default.odc) rather than overwriting each time. (Fixes a
+ bug that occurred when R was interrupted in the middle of a Bugs run.)
+ 2. Display different colored dots in the right panel of the graphical
+ display, to show the medians from each chain.
+Update 13 Mar 2003: fix minor bug in monitor()
+Update 10 Mar 2003: fix bug in pD and DIC calculations
+Update 7 Mar 2003:
+ 1. Fix display.parallel=T option by adding min.width so that very
+ intervals are still visible.
+ 2. Compute pD separately for each sequence (which gives much more
+ reasonable estimates before convergence).
+Update 8 Feb 2003: minor fixes in graphical display
+Update 6 Feb 2003:
+ 1. Approximate "effective sample size" n.eff given for each parameter.
+ 2. More explanatory material displayed.
+ 3. Use bringToTop() to automatically bring up the graphics window.
+Update 4 Feb 2003:
+ 1. Automatically compute the deviance, DIC, and pD. Bugs will not
+ always compute DIC and pD, so we do so using the definition,
+ DIC = E(deviance) + pD, using var(deviance)/2 as an estimate of pD.
+ (This is derived from the chi^2 distribution. We can't use the
+ Spiegelhalter et al. definition of DIC because we don't have access
+ to the deviance function.)
+ 2. Set default for n.thin so that, after thinning, the total number
+ saved iterations, n.sims, is approximately 1000.
+Update 14 Jan 2003 to run with the new WinBugs1.4. You may see an error
+ message and need to fix the dos.location assignment in bugs().
+Update 6 Jan 2003:
+ 1. Fix of bug that occurred with uppercase and lowercase variable names
+ 2. Set default for n.thin so that no more than about 500 iterations
+ will be saved from each sequence
+ 3. New option "display.parallel" added to show 80% inferences from
+ parallel sequences on the right panel of the graphical display. This
+ can be useful to understand what is going on when there are
+ convergence problems.
+Update 26 Dec 2002: fix of minwidth in bugs.plot.summary
+Update 11 Dec 2002:
+ 1. Automatic fixing of adaptive phases. Now you no longer need to run
+ for thousands of iterations when slice or Metropolis sampling is used.
+ 2. Various minor fixes
+Update 10 Dec 2002:
+ 1. Cool graphical display of convergence and inferences!
+ 2. New "attach2" function that overwrites so that all components of
+ the list are attached
+Update 29 Nov 2002:
+ 1. Fix of bug in 24 Nov update.
+ 2. Fix of bug in 16 Nov update.
+ 3. Length of chains is now pecified in terms of "n.iter" rather than
+ "n.keep".
+Update 24 Nov 2002: improved treatment of "parameters.to.save". For
+ example, you can now use "alpha" to indicate an entire array of parameters,
+ whereas before you had to save "alpha[]" or "alpha[,]" or whatever.
+Update 16 Nov 2002: mean, sd, median added to outputs
+Update 4 Nov 2002: more error-flagging added
+Update 26 Oct 2002:
+ 1. Parameters saved in order of the "parameters.to.save" vector
+ (not alphabetical order).
+ 2. Output saved in both matrix and list form.
+ 3. With the attach.sims=T setting (which is the default), the simulations
+ for all the saved parameters are saved as R objects. This is
+ convenient for later use of the simulations.
+Updates to 16 Oct 2002: more error-flagging added, mean/sd added to summary,
+ fixing scientific notation so Bugs can always read data and inits
+Update 21 Sept 2002: "quit=F" option changed to "debug=T"
+First version written 18 Sept 2002 by Andrew Gelman,
+ adapted from the EmBedBugs package by Kenneth Rice
Added: trunk/R2WinBUGS/DESCRIPTION
===================================================================
--- trunk/R2WinBUGS/DESCRIPTION (rev 0)
+++ trunk/R2WinBUGS/DESCRIPTION 2006-11-24 13:24:51 UTC (rev 3)
@@ -0,0 +1,19 @@
+Package: R2WinBUGS
+Title: Running WinBUGS and OpenBUGS from R
+Date: 2006-11-01
+Version: 2.0-4
+Author: originally written by Andrew Gelman <ge...@st...>;
+ changes and packaged by Sibylle Sturtz <st...@st...>
+ and Uwe Ligges <li...@st...>.
+ With considerable contributions by Gregor Gorjanc <gre...@bf...>
+ and Jouni Kerman <ke...@st...>.
+Description: Using this package,
+ it is possible to call a BUGS model, summarize inferences and
+ convergence in a table and graph, and save the simulations in arrays for easy access in R.
+Depends: R (>= 2.2.0)
+Suggests: coda (>= 0.9-0), BRugs (>= 0.3-0)
+SystemRequirements: WinBUGS 1.4 on Windows
+URL: http://www.stat.columbia.edu/~gelman/bugsR/
+Maintainer: Sibylle Sturtz <st...@st...>
+License: GPL version 2
+Packaged: Thu Nov 2 08:53:20 2006; ligges
Added: trunk/R2WinBUGS/NAMESPACE
===================================================================
--- trunk/R2WinBUGS/NAMESPACE (rev 0)
+++ trunk/R2WinBUGS/NAMESPACE 2006-11-24 13:24:51 UTC (rev 3)
@@ -0,0 +1,4 @@
+export(bugs, attach.all, detach.all, attach.bugs, detach.bugs, bugs.data, read.bugs,
+bugs.log, openbugs, as.bugs.array, write.model)
+S3method(print, bugs)
+S3method(plot, bugs)
Added: trunk/R2WinBUGS/R/as.bugs.array.R
===================================================================
--- trunk/R2WinBUGS/R/as.bugs.array.R (rev 0)
+++ trunk/R2WinBUGS/R/as.bugs.array.R 2006-11-24 13:24:51 UTC (rev 3)
@@ -0,0 +1,135 @@
+as.bugs.array <- function(sims.array, model.file=NULL, program=NULL, DIC=FALSE, n.iter=NULL, n.burnin=0, n.thin=1)
+{
+ # Jouni Kerman's function to convert a 3-way array to a Bugs object
+ #
+ # 'sims.array' is supposed to be a 3-way array with
+ # n.sims*n.chains*n.parameters simulations, and
+ # the 3rd component of dimnames(x) should have the parameter names.
+ d <- dim(sims.array)
+ n.keep <- d[1]
+ n.chains <- d[2]
+ n.parameters <- d[3]
+ n.sims <- n.keep*n.chains
+ if (is.null(n.iter)) n.iter <- (n.keep + n.burnin)*n.thin
+ #
+ parameter.names <- dimnames(sims.array)[[3]]
+ if (is.null(parameter.names)) {
+ parameter.names <- paste("P", 1:n.parameters, sep="")
+ dimnames(sims.array)[[3]] <- parameter.names
+ }
+ parameters.to.save <- unique(sapply(strsplit(parameter.names,"\\["),"[",1))
+ #
+ sims <- matrix(NA, n.sims, n.parameters)
+ root.long <- character(n.parameters)
+ indexes.long <- vector(n.parameters, mode = "list")
+ for (i in 1:n.parameters) {
+ temp <- decode.parameter.name(parameter.names[i])
+ root.long[i] <- temp$root
+ indexes.long[[i]] <- temp$indexes
+ }
+ n.roots <- length(parameters.to.save)
+ left.bracket.short <- as.vector(regexpr("[[]", parameters.to.save))
+ right.bracket.short <- as.vector(regexpr("[]]", parameters.to.save))
+ root.short <- ifelse(left.bracket.short == -1, parameters.to.save,
+ substring(parameters.to.save, 1, left.bracket.short - 1))
+ dimension.short <- rep(0, n.roots)
+ indexes.short <- vector(n.roots, mode = "list")
+ n.indexes.short <- vector(n.roots, mode = "list")
+ long.short <- vector(n.roots, mode = "list")
+ length.short <- numeric(n.roots)
+ for (j in 1:n.roots) {
+ long.short[[j]] <- (1:n.parameters)[root.long == root.short[j]]
+ length.short[j] <- length(long.short[[j]])
+ if (length.short[j] == 0){
+ stop(paste("parameter", root.short[[j]], "is not in the model"))
+ }
+ else if (length.short[j] > 1) {
+ dimension.short[j] <- length(indexes.long[[long.short[[j]][1]]])
+ n.indexes.short[[j]] <- numeric(dimension.short[j])
+ for (k in 1:dimension.short[j]){
+ n.indexes.short[[j]][k] <- length(unique(unlist(lapply(indexes.long[long.short[[j]]], .subset, k))))
+ }
+ length.short[j] <- prod(n.indexes.short[[j]])
+ if (length(long.short[[j]]) != length.short[j]){
+ stop(paste("error in parameter", root.short[[j]],
+ "in parameters.to.save"))
+ }
+ indexes.short[[j]] <- as.list(numeric(length.short[j]))
+ for (k in 1:length.short[j]){
+ indexes.short[[j]][[k]] <- indexes.long[[long.short[[j]][k]]]
+ }
+ }
+ }
+ rank.long <- unlist(long.short)
+ # -----
+ # yes, it's inefficient to do this, but for now I'm just letting this be as it is:
+ for (k in 1:n.parameters) {
+ sims[,k] <- as.vector(sims.array[,,k])
+ }
+ # ----
+ dimnames(sims) <- list(NULL, parameter.names)
+ summary <- monitor(sims.array, n.chains, keep.all = TRUE)
+ last.values <- as.list(numeric(n.chains))
+ for (i in 1:n.chains) {
+ n.roots.0 <- if (!is.null(DIC))
+ n.roots - 1
+ else n.roots
+ last.values[[i]] <- as.list(numeric(n.roots.0))
+ names(last.values[[i]]) <- root.short[1:n.roots.0]
+ for (j in 1:n.roots.0) {
+ if (dimension.short[j] <= 1) {
+ last.values[[i]][[j]] <- sims.array[n.keep, i, long.short[[j]]]
+ names(last.values[[i]][[j]]) <- NULL
+ }
+ else{
+ last.values[[i]][[j]] <- aperm(array(sims.array[n.keep, i,
+ long.short[[j]]], rev(n.indexes.short[[j]])), dimension.short[j]:1)
+ }
+ }
+ }
+ sims <- sims[sample(n.sims), ]
+ sims.list <- summary.mean <- summary.sd <- summary.median <- vector(n.roots,
+ mode = "list")
+ names(sims.list) <- names(summary.mean) <- names(summary.sd) <- names(summary.median) <- root.short
+ for (j in 1:n.roots) {
+ if (length.short[j] == 1) {
+ sims.list[[j]] <- sims[, long.short[[j]]]
+ summary.mean[[j]] <- summary[long.short[[j]], "mean"]
+ summary.sd[[j]] <- summary[long.short[[j]], "sd"]
+ summary.median[[j]] <- summary[long.short[[j]], "50%"]
+ }
+ else {
+ temp2 <- dimension.short[j]:1
+ sims.list[[j]] <- aperm(array(sims[, long.short[[j]]],
+ c(n.sims, rev(n.indexes.short[[j]]))), c(1, (dimension.short[j]+1):2))
+ summary.mean[[j]] <- aperm(array(summary[long.short[[j]],
+ "mean"], rev(n.indexes.short[[j]])), temp2)
+ summary.sd[[j]] <- aperm(array(summary[long.short[[j]],
+ "sd"], rev(n.indexes.short[[j]])), temp2)
+ summary.median[[j]] <- aperm(array(summary[long.short[[j]],
+ "50%"], rev(n.indexes.short[[j]])), temp2)
+ }
+ }
+ summary <- summary[rank.long, ]
+ all <- list(n.chains = n.chains, n.iter = n.iter, n.burnin = n.burnin,
+ n.thin = n.thin, n.keep = n.keep, n.sims = n.sims,
+ sims.array = sims.array[,,rank.long,drop = FALSE], sims.list = sims.list,
+ sims.matrix = sims[, rank.long], summary = summary, mean = summary.mean,
+ sd = summary.sd, median = summary.median, root.short = root.short,
+ long.short = long.short, dimension.short = dimension.short,
+ indexes.short = indexes.short, last.values = last.values, program=program,
+ model.file=model.file, is.DIC=!is.null(DIC), DIC=DIC)
+ if(sum(DIC)) {
+ deviance <- all$sims.array[, , dim(sims.array)[3], drop = FALSE]
+ dim(deviance) <- dim(deviance)[1:2]
+ pD <- numeric(n.chains)
+ DIC <- numeric(n.chains)
+ for (i in 1:n.chains) {
+ pD[i] <- var(deviance[, i])/2
+ DIC[i] <- mean(deviance[, i]) + pD[i]
+ }
+ all <- c(all, list(pD = mean(pD), DIC = mean(DIC)))
+ }
+ class(all) <- "bugs"
+ return(all)
+}
Property changes on: trunk/R2WinBUGS/R/as.bugs.array.R
___________________________________________________________________
Name: svn:keywords
+ Id
Added: trunk/R2WinBUGS/R/attach.all.R
===================================================================
--- trunk/R2WinBUGS/R/attach.all.R (rev 0)
+++ trunk/R2WinBUGS/R/attach.all.R 2006-11-24 13:24:51 UTC (rev 3)
@@ -0,0 +1,36 @@
+attach.all <- function(x, overwrite = NA, name = "attach.all"){
+ rem <- names(x) %in% ls(.GlobalEnv)
+ if(!any(rem)) overwrite <- FALSE
+ rem <- names(x)[rem]
+ if(is.na(overwrite)){
+ question <- paste("The following objects in .GlobalEnv will mask\nobjects in the attached database:\n",
+ paste(rem, collapse=", "),
+ "\nRemove these objects from .GlobalEnv?", sep="")
+ if(interactive()){
+ if(.Platform$OS.type == "windows")
+ overwrite <- "YES" == winDialog(type = "yesno", question)
+ else
+ overwrite <- 1 == menu(c("YES", "NO"), graphics = FALSE, title = question)
+ }
+ else overwrite <- FALSE
+ }
+ if(overwrite) remove(list=rem, envir=.GlobalEnv)
+ attach(x, name=name)
+}
+
+attach.bugs <- function (x, overwrite = NA){
+ if(class(x) != "bugs")
+ stop("attach.all() requires a bugs object.")
+ if("bugs.sims" %in% search()){
+ detach("bugs.sims")}
+ x$sims.list$n.sims <- x$n.sims # put n.sims into sims.list for convenience
+ r2 <- attach.all(x$sims.list, overwrite = overwrite, name = "bugs.sims")
+ invisible (bugs.sims = r2)
+}
+
+detach.all <- function(name = "attach.all")
+ do.call("detach", list(name=name))
+
+detach.bugs <- function(){
+ detach.all("bugs.sims")
+}
Property changes on: trunk/R2WinBUGS/R/attach.all.R
___________________________________________________________________
Name: svn:keywords
+ Id
Added: trunk/R2WinBUGS/R/bugs.R
===================================================================
--- trunk/R2WinBUGS/R/bugs.R (rev 0)
+++ trunk/R2WinBUGS/R/bugs.R 2006-11-24 13:24:51 UTC (rev 3)
@@ -0,0 +1,70 @@
+"bugs" <-
+function(data, inits, parameters.to.save, model.file = "model.bug",
+ n.chains = 3, n.iter = 2000, n.burnin = floor(n.iter / 2),
+ n.thin = max(1, floor(n.chains * (n.iter - n.burnin) / 1000)),
+ bin = (n.iter - n.burnin) / n.thin,
+ debug = FALSE, DIC = TRUE, digits = 5, codaPkg = FALSE,
+ bugs.directory = "c:/Program Files/WinBUGS14/",
+ program = c("winbugs", "openbugs", "WinBugs", "OpenBugs"),
+ working.directory = NULL,
+ clearWD = FALSE, useWINE = .Platform$OS.type != "windows", WINE = Sys.getenv("WINE"),
+ newWINE = FALSE, WINEPATH = NULL){
+
+ ## If OpenBUGS, we only call openbugs() and exit...
+ program <- match.arg(program)
+ if (program %in% c("openbugs", "OpenBugs"))
+ return(openbugs(data, inits, parameters.to.save, model.file,
+ n.chains, n.iter, n.burnin, n.thin, DIC, bugs.directory,
+ working.directory, digits))
+
+ ## Checking number of inits, which is NOT save here:
+ if(!missing(inits) && !is.function(inits) && !is.null(inits) && (length(inits) != n.chains))
+ stop("Number of initialized chains (length(inits)) != n.chains")
+ if (useWINE) { # attempt to find wine and winepath
+ if (!nchar(WINE)) {
+ WINE <- system("locate wine | grep bin/wine$", intern = TRUE)
+ WINE <- WINE[length(WINE)]
+ }
+ if (!length(WINE)) stop("couldn't locate WINE binary")
+ if (!nchar(WINEPATH)) {
+ WINEPATH <- system("locate winepath | grep bin/winepath$", intern = TRUE)
+ WINEPATH <- WINEPATH[length(WINEPATH)]
+ }
+ if (!length(WINEPATH)) stop("couldn't locate WINEPATH binary")
+ }
+ if(!is.null(working.directory)){
+ savedWD <- getwd()
+ setwd(working.directory)
+ on.exit(setwd(savedWD))
+ }
+ if(!file.exists(model.file)) stop(paste(model.file, "does not exist."))
+ if(file.info(model.file)$isdir) stop(paste(model.file, "is a directory, but a file is required."))
+ if(!(length(data) == 1 && is.vector(data) && is.character(data) && data == "data.txt")){
+ bugs.data(data, dir = getwd(), digits)
+ } else if(!file.exists(data))
+ stop("File data.txt does not exist.")
+ bugs.inits(inits, n.chains, digits)
+ if(DIC) parameters.to.save <- c(parameters.to.save, "deviance")
+ ## Model files with extension ".bug" need to be renamed to ".txt"
+ if(length(grep("[.]bug$", model.file))){
+ new.model.file <- sub("[.]bug$", ".txt", model.file)
+ file.copy(model.file, new.model.file, overwrite = TRUE)
+ on.exit(try(file.remove(new.model.file)), add = TRUE)
+ } else new.model.file <- model.file
+ bugs.script(parameters.to.save, n.chains, n.iter, n.burnin, n.thin,
+ new.model.file, debug=debug, is.inits=!is.null(inits),
+ bin = bin, DIC = DIC, useWINE = useWINE, newWINE = newWINE, WINEPATH = WINEPATH)
+ bugs.run(n.burnin, bugs.directory, WINE = WINE, useWINE = useWINE)
+
+ if(codaPkg)
+ return(file.path(getwd(), paste("coda", 1:n.chains, ".txt", sep="")))
+
+ sims <- c(bugs.sims(parameters.to.save, n.chains, n.iter, n.burnin, n.thin, DIC),
+ model.file = model.file, is.DIC = !is.null(DIC), program = program)
+ if(clearWD)
+ file.remove(c("data.txt", "log.odc", "log.txt", "codaIndex.txt",
+ paste("inits", 1:n.chains, ".txt", sep=""),
+ paste("coda", 1:n.chains, ".txt", sep="")))
+ class(sims) <- "bugs"
+ return(sims)
+}
Property changes on: trunk/R2WinBUGS/R/bugs.R
___________________________________________________________________
Name: svn:keywords
+ Id
Added: trunk/R2WinBUGS/R/bugs.data.R
===================================================================
--- trunk/R2WinBUGS/R/bugs.data.R (rev 0)
+++ trunk/R2WinBUGS/R/bugs.data.R 2006-11-24 13:24:51 UTC (rev 3)
@@ -0,0 +1,12 @@
+"bugs.data" <-
+function(data, dir = getwd(), digits = 5){
+ if(is.numeric(unlist(data)))
+ write.datafile(lapply(data, formatC, digits = digits, format = "E"),
+ file.path(dir, "data.txt"))
+ else {
+ data.list <- lapply(as.list(data), get, pos = parent.frame(2))
+ names(data.list) <- as.list(data)
+ write.datafile(lapply(data.list, formatC, digits = digits, format = "E"),
+ file.path(dir, "data.txt"))
+ }
+}
Property changes on: trunk/R2WinBUGS/R/bugs.data.R
___________________________________________________________________
Name: svn:keywords
+ Id
Added: trunk/R2WinBUGS/R/bugs.inits.R
===================================================================
--- trunk/R2WinBUGS/R/bugs.inits.R (rev 0)
+++ trunk/R2WinBUGS/R/bugs.inits.R 2006-11-24 13:24:51 UTC (rev 3)
@@ -0,0 +1,13 @@
+"bugs.inits" <-
+function (inits, n.chains, digits){
+ if(!is.null(inits)){
+ for (i in 1:n.chains){
+ if (is.function(inits))
+ write.datafile(lapply(inits(), formatC, digits = digits, format = "E"),
+ paste ("inits", i, ".txt", sep=""))
+ else
+ write.datafile(lapply(inits[[i]], formatC, digits = digits, format = "E"),
+ paste ("inits", i, ".txt", sep=""))
+ }
+ }
+}
Property changes on: trunk/R2WinBUGS/R/bugs.inits.R
___________________________________________________________________
Name: svn:keywords
+ Id
Added: trunk/R2WinBUGS/R/bugs.log.R
===================================================================
--- trunk/R2WinBUGS/R/bugs.log.R (rev 0)
+++ trunk/R2WinBUGS/R/bugs.log.R 2006-11-24 13:24:51 UTC (rev 3)
@@ -0,0 +1,46 @@
+bugs.log <- function(file)
+{
+ if(!file.exists(file))
+ stop("Log file", file, "does not exist")
+ logfile <- readLines(file)
+
+ statsStart <- which(logfile == "Node statistics") + 2
+ if(!length(statsStart))
+ stop("Log file", file, "does not contain node statistics.")
+ ## + 2 to remove
+ ## "Node statistics"
+ ## "\t node\t mean\t sd\t MC error\t2.5%\tmedian\t97.5%\tstart\tsample"
+
+ statsEnd <- which(logfile == "dic.stats()") - 1
+ ## + 1 to remove
+ ## "dic.stats()"
+
+ statsTable <- logfile[statsStart:statsEnd]
+ statsTable <- sub("\t", "", statsTable)
+ statsTable <- sapply(strsplit(statsTable, "\t"), "[")
+ colnames(statsTable) <- statsTable[1,]
+ statsTable <- t(apply(statsTable[-1,], 2, as.numeric))
+ colnames(statsTable) <- c("mean", "sd", "MC error", "2.5%", "median", "97.5%", "start", "sample")
+
+ ## DIC
+ DICStart <- which(logfile == "DIC") + 3
+ ## + 3 to remove
+ ## "DIC"
+ ## "Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes"
+ ## "\tDbar\tDhat\tpD\tDIC\t"
+
+ DICEnd <- grep("history(", logfile, fixed = TRUE) - 1
+ ## - 1 to remove
+ ## "history(+, ..."
+
+ if(!length(DICEnd) || !length(DICStart) || (DICEnd < DICStart)){
+ DICTable <- NA
+ }else{
+ DICTable <- logfile[DICStart:DICEnd]
+ DICTable <- sapply(strsplit(DICTable, "\t"), "[")
+ colnames(DICTable) <- DICTable[1,]
+ DICTable <- t(apply(DICTable[-1,], 2, as.numeric))
+ colnames(DICTable) <- c("Dbar", "Dhat", "pD", "DIC")
+ }
+ return(list(stats = statsTable, DIC = DICTable))
+}
Property changes on: trunk/R2WinBUGS/R/bugs.log.R
___________________________________________________________________
Name: svn:keywords
+ Id
Added: trunk/R2WinBUGS/R/bugs.plot.inferences.R
===================================================================
--- trunk/R2WinBUGS/R/bugs.plot.inferences.R (rev 0)
+++ trunk/R2WinBUGS/R/bugs.plot.inferences.R 2006-11-24 13:24:51 UTC (rev 3)
@@ -0,0 +1,110 @@
+"bugs.plot.inferences" <-
+function (sims, display.parallel, ...){
+ if (.Device=="windows" ||
+ (.Device=="null device" && options("device")=="windows")){
+ cex.names <- .7
+ cex.axis <- .6
+ cex.tiny <- .4
+ cex.points <- .7
+ standard.width <- 30
+ max.width <- 40
+ min.width <- .02
+ }
+ else {
+ cex.names <- .7
+ cex.axis <- .6
+ cex.tiny <- .4
+ cex.points <- .3
+ standard.width <- 30
+ max.width <- 40
+ min.width <- .01
+ }
+ rootnames <- sims$root.short
+ n.roots <- length(rootnames)
+ sims.array <- sims$sims.array
+ n.chains <- sims$n.chains
+ dimension.short <- sims$dimension.short
+ indexes.short <- sims$indexes.short
+ long.short <- sims$long.short
+ height <- .6
+ par (mar=c(0,0,1,0))
+ plot (c(0,1), c(-n.roots-.5,-.4),
+ ann=FALSE, bty="n", xaxt="n", yaxt="n", type="n")
+ W <- max(strwidth(rootnames, cex=cex.names))
+ B <- (1-W)/3.8
+ A <- 1-3.5*B
+ if (display.parallel)
+ text (A, -.4, "80% interval for each chain", adj=0, cex=cex.names)
+ else
+ text (A, -.4, "medians and 80% intervals", adj=0, cex=cex.names)
+ num.height <- strheight (1:9, cex=cex.tiny)
+ for (k in 1:n.roots){
+ text (0, -k, rootnames[k], adj=0, cex=cex.names)
+ J <- min (length(long.short[[k]]), max.width)
+ if (k==1)
+ index <- 1:J
+ else
+ index <- sum (unlist(lapply(long.short,length))[1:(k-1)]) + 1:J
+ spacing <- 3.5/max(J,standard.width)
+ med <- numeric(J)
+ i80 <- matrix( , J, 2)
+ med.chains <- matrix( , J, sims$n.chains)
+ i80.chains <- array(NA, c(J, sims$n.chains, 2))
+ for (j in 1:J){
+ med[j] <- median(sims.array[,,index[j]])
+ i80[j,] <- quantile(sims.array[,,index[j]], c(.1,.9))
+ for (m in 1:n.chains){
+ med.chains[j,m] <- quantile (sims.array[,m,index[j]], .5)
+ i80.chains[j,m,] <- quantile (sims.array[,m,index[j]], c(.1,.9))
+ }
+ }
+ rng <- range (i80, i80.chains)
+ p.rng <- pretty(rng, n = 2)
+ b <- height/(max(p.rng) - min(p.rng))
+ a <- -(k+height/2) - b*p.rng[1]
+ lines (A+c(0,0), -k+c(-height/2,height/2))
+#
+# plot a line at zero (if zero is in the range of the mini-plot)
+#
+ if (min(p.rng)<0 & max(p.rng)>0)
+ lines (A+B*spacing*c(0,J+1), rep (a,2), lwd=.5, col="gray")
+
+ for (x in p.rng){
+ text (A-B*.2, a+b*x, x, cex=cex.axis)
+ lines (A+B*c(-.05,0), rep(a+b*x,2))
+ }
+ for (j in 1:J){
+ if (display.parallel){
+ for (m in 1:n.chains){
+ interval <- a + b*i80.chains[j,m,]
+ if (interval[2]-interval[1] < min.width)
+ interval <- mean(interval) + c(-1,1)*min.width/2
+ lines (A+B*spacing*rep(j+.6*(m-(n.chains+1)/2)/n.chains,2),
+ interval, lwd=.5, col=m+1)
+ }
+ }
+ else {
+ lines (A+B*spacing*rep(j,2), a + b*i80[j,], lwd=.5)
+ for (m in 1:n.chains)
+# points (A+B*spacing*j, a + b*med[j], pch=20, cex=cex.points)
+ points (A+B*spacing*j, a + b*med.chains[j,m], pch=20, cex=cex.points,
+ col=m+1)
+ }
+ dk <- dimension.short[k]
+ if (dk>0){
+ for (m in 1:dk){
+ index0 <- indexes.short[[k]][[j]][m]
+ if (j==1)
+ text(A+B*spacing*j, -k-height/2-.05-num.height*(m-1), index0,
+ cex=cex.tiny)
+ else if (index0!=indexes.short[[k]][[j-1]][m] &
+ (index0%%(floor(log10(index0)+1))==0))
+ text(A+B*spacing*j, -k-height/2-.05-num.height*(m-1), index0,
+ cex=cex.tiny)
+ }
+ }
+ }
+ if (J<length(long.short[[k]])) text (-.015, -k, "*",
+ cex=cex.names, col="red")
+ }
+}
Property changes on: trunk/R2WinBUGS/R/bugs.plot.inferences.R
___________________________________________________________________
Name: svn:keywords
+ Id
Added: trunk/R2WinBUGS/R/bugs.plot.summary.R
===================================================================
--- trunk/R2WinBUGS/R/bugs.plot.summary.R (rev 0)
+++ trunk/R2WinBUGS/R/bugs.plot.summary.R 2006-11-24 13:24:51 UTC (rev 3)
@@ -0,0 +1,124 @@
+"bugs.plot.summary" <-
+function (sims, ...){
+ DIC <- sims$is.DIC
+ if (.Device=="windows" ||
+ (.Device=="null device" && options("device")=="windows")){
+ cex.names <- .7
+ cex.top <- .7
+ cex.points <- .7
+ max.length <- 50
+ min.width <- .01
+ }
+ else {
+ cex.names <- .7
+ cex.top <- .7
+ cex.points <- .3
+ max.length <- 80
+ min.width <- .005
+ }
+ summ <- sims$summary
+ sims.array <- sims$sims.array
+ n.chains <- sims$n.chains
+ n.parameters <- nrow(summ)
+
+ J0 <- unlist(lapply(sims$long.short, length))
+ if (DIC) J0 <- J0[1:(length(J0)-1)] # don't display deviance summaries
+ J <- J0
+ total <- ceiling(sum(J+.5))
+ while ((total > max.length) && max(J)>1){### vielleicht optimieren ...
+ J[J==max(J)] <- max(J)-1
+ total <- ceiling(sum(J+.5))
+ }
+ pos <- -1
+ ypos <- NULL
+ id <- NULL
+ ystart <- NULL
+ jj <- 1:J[1]
+ n.roots <- length(sims$root.short)
+ if (DIC) n.roots <- n.roots-1 # don't display deviance summaries
+ ystart <- numeric(n.roots)
+ for (k in 1:n.roots){
+ ystart[k] <- pos
+ ypos <- c(ypos, pos - seq(0, J[k]-1))
+ id <- c(id, 1:J[k])
+ pos <- pos - J[k] -.5
+ if (k>1) jj <- c(jj, sum(J0[1:(k-1)]) + (1:J[k]))
+ }
+ bottom <- min(ypos)-1
+
+ med <- numeric(sum(J))
+ i80 <- matrix( , sum(J), 2)
+ i80.chains <- array (NA, c(sum(J), n.chains, 2))
+ for (j in 1:sum(J)){
+ med[j] <- median (sims.array[,,jj[j]])
+ i80[j,] <- quantile (sims.array[,,jj[j]], c(.1,.9))
+ for (m in 1:n.chains)
+ i80.chains[j,m,] <- quantile (sims.array[,m,jj[j]], c(.1,.9))
+ }
+ rng <- range (i80, i80.chains)
+ p.rng <- pretty(rng, n = 2)
+ b <- 2 / (max(p.rng) - min(p.rng))
+ a <- -b * p.rng[1]
+
+ par (mar=c(0,0,1,3))
+ plot (c(0,1), c(min(bottom, -max.length)-3,2.5),
+ ann=FALSE, bty="n", xaxt="n", yaxt="n", type="n")
+ W <- max(strwidth(unlist(dimnames(summ)[[1]]), cex=cex.names))
+ B <- (1-W)/3.6
+ A <- 1-3.5*B
+ B <- (1-A)/3.5
+ b <- B*b
+ a <- A + B*a
+ text (A+B*1, 2.5, "80% interval for each chain", cex=cex.top)
+ lines (A+B*c(0,2), c(0,0))
+ lines (A+B*c(0,2), rep(bottom,2))
+ if(n.chains > 1){
+ text (A+B*3, 2.6, "R-hat", cex=cex.top)
+ lines (A+B*c(2.5,3.5), c(0,0))
+ lines (A+B*c(2.5,3.5), rep(bottom,2))
+ }
+#
+# line at zero
+#
+ if (min(p.rng)<0 & max(p.rng)>0)
+ lines (rep(a,2), c(0,bottom), lwd=.5, col="gray")
+
+ for (x in p.rng){
+ text (a+b*x, 1, x, cex=cex.names)
+ lines (rep(a+b*x,2), c(0,-.2))
+ text (a+b*x, bottom-1, x, cex=cex.names)
+ lines (rep(a+b*x,2), bottom+c(0,.2))
+ }
+ if(n.chains > 1)
+ for (x in seq(1,2,.5)){
+ text (A+B*(1.5+seq(1,2,.5)), rep(1,3), c("1","1.5","2+"), cex=cex.names)
+ lines (A+B*rep(1.5+x,2), c(0,-.2))
+ text (A+B*(1.5+seq(1,2,.5)), rep(bottom-1,3), c("1","1.5","2+"),
+ cex=cex.names)
+ lines (A+B*rep(1.5+x,2), bottom+c(0,.2))
+ }
+ for (j in 1:sum(J)){
+ name <- dimnames(summ)[[1]][jj[j]]
+ if (id[j]==1)
+ text (0, ypos[j], name, adj=0, cex=cex.names)
+ else {
+ pos <- as.vector(regexpr("[[]", name))
+ text (strwidth(substring(name,1,pos-1),cex=cex.names),
+ ypos[j], substring(name, pos, nchar(name)), adj=0, cex=cex.names)
+ }
+ for (m in 1:n.chains){
+ interval <- a + b*i80.chains[j,m,]
+ if (interval[2]-interval[1] < min.width)
+ interval <- mean(interval) + c(-1,1)*min.width/2
+ lines (interval, rep(ypos[j]-.1*(m-(n.chains+1)/2),2), lwd=1, col=m+1)
+ if(n.chains > 1)
+ points (A+B*(1.5 + min(max(summ[jj[j],"Rhat"],1),2)), ypos[j], pch=20, cex=cex.points)
+ }
+ }
+ for (k in 1:n.roots){
+ if (J[k]<J0[k]) text (-.015, ystart[k], "*", cex=cex.names,
+ col="red")
+ }
+ if (sum(J!=J0)>0) text (0, bottom-3,
+ "* array truncated for lack of space", adj=0, cex=cex.names, col="red")
+}
Property changes on: trunk/R2WinBUGS/R/bugs.plot.summary.R
___________________________________________________________________
Name: svn:keywords
+ Id
Added: trunk/R2WinBUGS/R/bugs.run.R
===================================================================
--- trunk/R2WinBUGS/R/bugs.run.R (rev 0)
+++ trunk/R2WinBUGS/R/bugs.run.R 2006-11-24 13:24:51 UTC (rev 3)
@@ -0,0 +1,35 @@
+"bugs.run" <-
+ function(n.burnin, bugs.directory, WINE = "",
+ useWINE = .Platform$OS.type != "windows", newWINE = TRUE){
+
+if(useWINE && !newWINE) bugs.directory <- win2native(bugs.directory)
+
+## Update the lengths of the adaptive phases in the Bugs updaters
+ try(bugs.update.settings(n.burnin, bugs.directory))
+## Return the lengths of the adaptive phases to their original settings
+ on.exit(try(file.copy(file.path(bugs.directory, "System/Rsrc/Registry_Rsave.odc"),
+ file.path(bugs.directory, "System/Rsrc/Registry.odc"),
+ overwrite = TRUE)))
+## Search Win*.exe (WinBUGS executable) within bugs.directory
+ dos.location <- file.path(bugs.directory,
+ grep("^Win[[:alnum:]]*[.]exe$", list.files(bugs.directory), value = TRUE)[1])
+ if(!file.exists(dos.location))
+ stop(paste("WinBUGS executable does not exist in", bugs.directory))
+## Call Bugs and have it run with script.txt
+ bugsCall <- paste("\"", dos.location, "\" /par \"",
+ native2win(file.path(getwd(), "script.txt"), newWINE = newWINE),
+ "\"", sep = "")
+ if (useWINE)
+ bugsCall <- paste(WINE, bugsCall)
+ temp <- system(bugsCall)
+
+ if(temp == -1)
+ stop("Error in bugs.run().\nCheck that WinBUGS is in the specified directory.")
+## Stop and print an error message if Bugs did not run correctly
+ if (length(grep("Bugs did not run correctly",
+ scan("coda1.txt", character(), quiet=TRUE, sep="\n"))) > 0)
+ stop("Look at the log file and\ntry again with debug=TRUE and figure out what went wrong within Bugs.")
+}
+
+
+
Property changes on: trunk/R2WinBUGS/R/bugs.run.R
___________________________________________________________________
Name: svn:keywords
+ Id
Added: trunk/R2WinBUGS/R/bugs.script.R
===================================================================
--- trunk/R2WinBUGS/R/bugs.script.R (rev 0)
+++ trunk/R2WinBUGS/R/bugs.script.R 2006-11-24 13:24:51 UTC (rev 3)
@@ -0,0 +1,48 @@
+"bugs.script" <-
+function (parameters.to.save, n.chains, n.iter, n.burnin,
+ n.thin, model.file, debug=FALSE, is.inits, bin,
+ DIC = FALSE, useWINE = FALSE, newWINE = FALSE, WINEPATH = NULL){
+### Write file script.txt for Bugs to read
+### if (n.chains<2) stop ("n.chains must be at least 2")
+
+ if((ceiling(n.iter/n.thin) - ceiling(n.burnin/n.thin)) < 2)
+ stop ("(n.iter-n.burnin)/n.thin must be at least 2")
+ working.directory <- getwd()
+ script <- "script.txt"
+ model <- if(length(grep("\\\\", model.file)) || length(grep("/", model.file))){
+ model.file
+ } else file.path(working.directory, model.file)
+ data <- file.path(working.directory, "data.txt")
+ history <- file.path(working.directory, "history.odc")
+ coda <- file.path(working.directory, "coda")
+ logfile <- file.path(working.directory, "log.odc")
+ logfileTxt <- file.path(working.directory, "log.txt")
+ inits <- sapply(paste(working.directory, "/inits", 1:n.chains, ".txt", sep=""), native2win)
+ initlist <- paste("inits (", 1:n.chains, ", '", inits, "')\n", sep="")
+ savelist <- paste("set (", parameters.to.save, ")\n", sep="")
+ redo <- ceiling((n.iter-n.burnin)/(n.thin*bin))
+ cat(
+ "display ('log')\n",
+ "check ('", native2win(model), "')\n",
+ "data ('", native2win(data), "')\n",
+ "compile (", n.chains, ")\n",
+ if(is.inits) initlist,
+ "gen.inits()\n",
+ "thin.updater (", n.thin, ")\n",
+ "update (", ceiling(n.burnin/n.thin), ")\n",
+ savelist,
+ if(DIC) "dic.set()\n",
+ rep(
+ c("update (", formatC(ceiling(bin), format = "d"), ")\n",
+ "coda (*, '", native2win(coda), "')\n"),redo),
+ "stats (*)\n",
+ if(DIC) "dic.stats()\n",
+ "history (*, '", native2win(history), "')\n",
+ "save ('", native2win(logfile), "')\n",
+ "save ('", native2win(logfileTxt), "')\n",
+ file=script, sep="", append=FALSE)
+ if (!debug) cat ("quit ()\n", file=script, append=TRUE)
+ sims.files <- paste ("coda", 1:n.chains, ".txt", sep="")
+ for (i in 1:n.chains) cat ("WinBUGS did not run correctly.\n",
+ file=sims.files[i], append=FALSE)
+}
Property changes on: trunk/R2WinBUGS/R/bugs.script.R
___________________________________________________________________
Name: svn:keywords
+ Id
Added: trunk/R2WinBUGS/R/bugs.sims.R
===================================================================
--- trunk/R2WinBUGS/R/bugs.sims.R (rev 0)
+++ trunk/R2WinBUGS/R/bugs.sims.R 2006-11-24 13:24:51 UTC (rev 3)
@@ -0,0 +1,130 @@
+"bugs.sims" <-
+function (parameters.to.save, n.chains, n.iter, n.burnin, n.thin, DIC = TRUE){
+# Read the simulations from Bugs into R, format them, and monitor convergence
+ sims.files <- paste ("coda", 1:n.chains, ".txt", sep="")
+ index <- read.table ("codaIndex.txt", header=FALSE, sep="\t")
+ parameter.names <- as.vector(index[,1])
+ n.keep <- index[1,3] - index[1,2] + 1
+ n.parameters <- length(parameter.names)
+ n.sims <- n.keep*n.chains
+ sims <- matrix( , n.sims, n.parameters)
+ sims.array <- array (NA, c(n.keep, n.chains, n.parameters))
+ root.long <- character(n.parameters)
+ indexes.long <- vector(n.parameters, mode = "list")
+ for (i in 1:n.parameters){
+ temp <- decode.parameter.name(parameter.names[i])
+ root.long[i] <- temp$root
+ indexes.long[[i]] <- temp$indexes
+ }
+ n.roots <- length(parameters.to.save)
+ left.bracket.short <- as.vector (regexpr("[[]", parameters.to.save))
+ right.bracket.short <- as.vector (regexpr("[]]", parameters.to.save))
+ root.short <- ifelse (left.bracket.short==-1, parameters.to.save,
+ substring (parameters.to.save, 1, left.bracket.short-1))
+ dimension.short <- rep(0, n.roots)
+ indexes.short <- vector(n.roots, mode = "list")
+ n.indexes.short <- vector(n.roots, mode = "list")
+ long.short <- vector(n.roots, mode = "list")
+ length.short <- numeric(n.roots)
+ ##SS, UL##: Let's optimize the following loops ...
+ for (j in 1:n.roots){
+ long.short[[j]] <- (1:n.parameters)[root.long==root.short[j]]
+ length.short[j] <- length(long.short[[j]])
+ if (length.short[j]==0)
+ stop (paste ("parameter", root.short[[j]], "is not in the model"))
+ else if (length.short[j]>1){
+ dimension.short[j] <- length(indexes.long[[long.short[[j]][1]]])
+ n.indexes.short[[j]] <- numeric(dimension.short[j])
+ for (k in 1:dimension.short[j]) n.indexes.short[[j]][k] <- length (
+ unique (unlist (lapply (indexes.long[long.short[[j]]], .subset, k))))
+ length.short[j] <- prod(n.indexes.short[[j]])
+ if (length(long.short[[j]])!=length.short[j]) stop (paste
+ ("error in parameter", root.short[[j]], "in parameters.to.save"))
+ indexes.short[[j]] <- as.list(numeric(length.short[j]))
+ for (k in 1:length.short[j])
+ indexes.short[[j]][[k]] <- indexes.long[[long.short[[j]][k]]]
+ }
+ }
+# rank.long <- rank(paste(rep(root.short,length.short),
+# (1:n.parameters)/10^ceiling(log10(n.parameters)),sep="."))
+ rank.long <- unlist(long.short)
+
+ for (i in 1:n.chains){
+ sims.i <- scan (sims.files[i], quiet=TRUE) [2*(1:(n.keep*n.parameters))]
+ sims[(n.keep*(i-1)+1):(n.keep*i), ] <- sims.i
+ sims.array[,i,] <- sims.i
+ }
+ dimnames (sims) <- list (NULL, parameter.names)
+ dimnames (sims.array) <- list (NULL, NULL, parameter.names)
+#
+# Perform convergence checks and compute medians and quantiles.
+#
+ summary <- monitor (sims.array, n.chains, keep.all=TRUE)
+#
+# Create outputs
+#
+ last.values <- as.list (numeric(n.chains))
+ for (i in 1:n.chains){
+ n.roots.0 <- if(DIC) n.roots-1 else n.roots
+ last.values[[i]] <- as.list (numeric(n.roots.0))
+ names(last.values[[i]]) <- root.short[1:n.roots.0]
+ for (j in 1:n.roots.0){
+ if (dimension.short[j]<=1){
+ last.values[[i]][[j]] <- sims.array[n.keep,i,long.short[[j]]]
+ names(last.values[[i]][[j]]) <- NULL
+ }
+ else
+ last.values[[i]][[j]] <- aperm (array(sims.array[n.keep,i,long.short[[j]]],
+ rev(n.indexes.short[[j]])), dimension.short[j]:1)
+ }
+ }
+ sims <- sims [sample(n.sims),] # scramble (for convenience in analysis)
+ sims.list <- summary.mean <- summary.sd <- summary.median <- vector(n.roots, mode = "list")
+ names(sims.list) <- names(summary.mean) <- names(summary.sd) <- names(summary.median) <- root.short
+ for (j in 1:n.roots){
+ if (length.short[j]==1){
+ sims.list[[j]] <- sims[,long.short[[j]]]
+ summary.mean[[j]] <- summary[long.short[[j]],"mean"]
+ summary.sd[[j]] <- summary[long.short[[j]],"sd"]
+ summary.median[[j]] <- summary[long.short[[j]],"50%"]
+ }
+ else{
+ temp2 <- dimension.short[j]:1
+ sims.list[[j]] <- aperm (array (sims[,long.short[[j]]],
+ c(n.sims,rev(n.indexes.short[[j]]))), c(1,(dimension.short[j]+1):2))
+ summary.mean[[j]] <- aperm (array (summary[long.short[[j]],"mean"],
+ rev(n.indexes.short[[j]])), temp2)
+ summary.sd[[j]] <- aperm (array (summary[long.short[[j]],"sd"],
+ rev(n.indexes.short[[j]])), temp2)
+ summary.median[[j]] <- aperm (array (summary[long.short[[j]],"50%"],
+ rev(n.indexes.short[[j]])), temp2)
+ }
+ }
+ summary <- summary[rank.long,]
+ all <- list (n.chains=n.chains, n.iter=n.iter, n.burnin=n.burnin,
+ n.thin=n.thin, n.keep=n.keep, n.sims=n.sims,
+ sims.array=sims.array[,,rank.long,drop=FALSE], sims.list=sims.list,
+ sims.matrix=sims[,rank.long], summary=summary, mean=summary.mean,
+ sd=summary.sd, median=summary.median, root.short=root.short,
+ long.short=long.short, dimension.short=dimension.short,
+ indexes.short=indexes.short, last.values=last.values)
+ if(DIC){
+ LOG <- bugs.log("log.txt")$DIC
+ if(any(is.na(LOG))){
+ deviance <- all$sims.array[, , dim(sims.array)[3], drop = FALSE]
+ dim(deviance) <- dim(deviance)[1:2]
+ pD <- numeric(n.chains)
+ DIC <- numeric(n.chains)
+ for (i in 1:n.chains) {
+ pD[i] <- var(deviance[, i])/2
+ DIC[i] <- mean(deviance[, i]) + pD[i]
+ }
+ all <- c(all, list(pD = mean(pD), DIC = mean(DIC), DICbyR = TRUE))
+ } else {
+ DIC <- LOG[nrow(LOG),4]
+ pD <- LOG[nrow(LOG),3]
+ all <- c(all, list(pD = pD, DIC = DIC, DICbyR = FALSE))
+ }
+ }
+ return(all)
+}
Property changes on: trunk/R2WinBUGS/R/bugs.sims.R
___________________________________________________________________
Name: svn:keywords
+ Id
Added: trunk/R2WinBUGS/R/bugs.update.settings.R
===================================================================
--- trunk/R2WinBUGS/R/bugs.update.settings.R (rev 0)
+++ trunk/R2WinBUGS/R/bugs.update.settings.R 2006-11-24 13:24:51 UTC (rev 3)
@@ -0,0 +1,30 @@
+"bugs.update.settings" <-
+function (n.burnin, bugs.directory){
+ char.burnin <- as.character(n.burnin - 1)
+ file.copy(file.path(bugs.directory, "System/Rsrc/Registry.odc"),
+ file.path(bugs.directory, "System/Rsrc/Registry_Rsave.odc"),
+ overwrite = TRUE)
+ registry <- readBin(file.path(bugs.directory, "System/Rsrc/Registry.odc"),
+ "character", 400, size = 1, endian = "little")
+ locale <- Sys.getlocale("LC_CTYPE")
+ Sys.setlocale("LC_CTYPE", "C")
+ info <- registry[regexpr("Int", registry, fixed = TRUE, useBytes = TRUE) > 0]
+ while(regexpr("\r", info) > 0){
+ newline <- regexpr("\r", info)
+ info <- substring(info, newline + 1)
+ line <- substring(info, 1, regexpr("\r", info) - 1)
+ if(regexpr("AdaptivePhase", line) > 0){
+ numpos <- regexpr("Int", line, fixed = TRUE, useBytes = TRUE) + 4
+ num <- substring(line, numpos)
+ if (as.numeric(num) > n.burnin){
+ blanks <- rep(" ", nchar(num, type = "chars") - nchar(char.burnin, type = "chars"))
+ num.new <- paste(paste(blanks, collapse = ""), char.burnin, sep = "")
+ line.new <- sub(num, num.new, line)
+ registry <- sub(line, line.new, registry)
+ }
+ }
+ }
+ Sys.setlocale("LC_CTYPE", locale)
+ writeBin (registry,
+ file.path(bugs.directory, "System/Rsrc/Registry.odc"), endian = "little")
+}
Property changes on: trunk/R2WinBUGS/R/bugs.update.settings.R
___________________________________________________________________
Name: svn:keywords
+ Id
Added: trunk/R2WinBUGS/R/conv.par.R
===================================================================
--- trunk/R2WinBUGS/R/conv.par.R (rev 0)
+++ trunk/R2WinBUGS/R/conv.par.R 2006-11-24 13:24:51 UTC (rev 3)
@@ -0,0 +1,86 @@
+"conv.par" <-
+function (x, n.chains, Rupper.keep = TRUE) {
+ m <- ncol(x)
+ n <- nrow(x)
+
+# We compute the following statistics:
+#
+# xdot: vector of sequence means
+# s2: vector of sequence sample variances (dividing by n-1)
+# W = mean(s2): within MS
+# B = n*var(xdot): between MS.
+# muhat = mean(xdot): grand mean; unbiased under strong stationarity
+# varW = var(s2)/m: estimated sampling var of W
+# varB = B^2 * 2/(m+1): estimated sampling var of B
+# covWB = (n/m)*(cov(s2,xdot^2) - 2*muhat*cov(s^2,xdot)):
+# estimated sampling cov(W,B)
+# sig2hat = ((n-1)/n))*W + (1/n)*B: estimate of sig2; unbiased under
+# strong stationarity
+# quantiles: emipirical quantiles from last half of simulated sequences
+
+ xdot <- apply(x,2,mean)
+ muhat <- mean(xdot)
+ s2 <- apply(x,2,var)
+ W <- mean(s2)
+ quantiles <- quantile (as.vector(x), probs=c(.025,.25,.5,.75,.975))
+
+ if ((W > 1.e-8) && (n.chains > 1)) { # non-degenerate case
+
+ B <- n*var(xdot)
+ varW <- var(s2)/m
+ varB <- B^2 * 2/(m-1)
+ covWB <- (n/m)*(cov(s2,xdot^2) - 2*muhat*cov(s2,xdot))
+ sig2hat <- ((n-1)*W + B)/n
+
+# Posterior interval post.range combines all uncertainties
+# in a t interval with center muhat, scale sqrt(postvar),
+# and postvar.df degrees of freedom.
+#
+# postvar = sig2hat + B/(mn): variance for the posterior interval
+# The B/(mn) term is there because of the
+# sampling variance of muhat.
+# varpostvar: estimated sampling variance of postvar
+
+ postvar <- sig2hat + B/(m*n)
+ varpostvar <- max(0,
+ (((n-1)^2) * varW + (1 + 1/m)^2 * varB + 2 * (n-1) * (1 + 1/m) * covWB) / n^2)
+ post.df <- min(2*(postvar^2/varpostvar), 1000)
+
+# Estimated potential scale reduction (that would be achieved by
+# continuing simulations forever) has two components: an estimate and
+# an approx. 97.5% upper bound.
+#
+# confshrink = sqrt(postvar/W),
+# multiplied by sqrt(df/(df-2)) as an adjustment for the
+### CHANGED TO sqrt((df+3)/(df+1))
+# width of the t-interval with df degrees of freedom.
+#
+# postvar/W = (n-1)/n + (1+1/m)(1/n)(B/W); we approximate the sampling dist.
+# of (B/W) by an F distribution, with degrees of freedom estimated
+# from the approximate chi-squared sampling dists for B and W. (The
+# F approximation assumes that the sampling dists of B and W are independent;
+# if they are positively correlated, the approximation is conservative.)
+
+ confshrink.range <- postvar/W
+ if(Rupper.keep){
+ varlo.df <- 2*(W^2/varW)
+ confshrink.range <- c(confshrink.range,
+ (n-1)/n + (1+1/m)*(1/n)*(B/W) * qf(.975, m-1, varlo.df))
+ }
+ confshrink.range <- sqrt(confshrink.range * (post.df+3) / (post.df+1))
+
+# Calculate effective sample size: m*n*min(sigma.hat^2/B,1)
+# This is a crude measure of sample size because it relies on the between
+# variance, B, which can only be estimated with m degrees of freedom.
+
+ n.eff <- m*n*min(sig2hat/B,1)
+ list(quantiles=quantiles, confshrink=confshrink.range,
+ n.eff=n.eff)
+
+ }
+ else { # degenerate case: all entries in "data matrix" are identical
+ list (quantiles=quantiles, confshrink = rep(1, Rupper.keep + 1),
+ n.eff=1)
+
+ }
+}
Property changes on: trunk/R2WinBUGS/R/conv.par.R
___________________________________________________________________
Name: svn:keywords
+ Id
Added: trunk/R2WinBUGS/R/decode.parameter.name.R
===================================================================
--- trunk/R2WinBUGS/R/decode.parameter.name.R (rev 0)
+++ trunk/R2WinBUGS/R/decode.parameter.name.R 2006-11-24 13:24:51 UTC (rev 3)
@@ -0,0 +1,22 @@
+"decode.parameter.name" <-
+function (a){
+#
+# Decodes Bugs parameter names
+# (e.g., "beta[3,14]" becomes "beta" with 2 indexes: 3 and 14)
+# for use by the bugs.sim() function
+#
+ left.bracket <- regexpr ("[[]", a)
+ if (left.bracket==-1){
+ root <- a
+ dimension <- 0
+ indexes <- NA
+ }
+ else {
+ root <- substring (a, 1, left.bracket-1)
+ right.bracket <- regexpr ("[]]", a)
+ a <- substring (a, left.bracket+1, right.bracket-1)
+ indexes <- as.numeric(unlist(strsplit(a, ",")))
+ dimension <- length(indexes)
+ }
+ return (list (root=root, dimension=dimension, indexes=indexes))
+}
Property changes on: trunk/R2WinBUGS/R/decode.parameter.name.R
___________________________________________________________________
Name: svn:keywords
+ Id
Added: trunk/R2WinBUGS/R/formatdata.R
===================================================================
--- trunk/R2WinBUGS/R/formatdata.R (rev 0)
+++ trunk/R2WinBUGS/R/formatdata.R 2006-11-24 13:24:51 UTC (rev 3)
@@ -0,0 +1,24 @@
+"formatdata" <-
+function (datalist){
+ if (!is.list(datalist) || is.data.frame(datalist))
+ stop("Argument to formatdata() must be a list.")
+ n <- length(datalist)
+ datalist.string <- vector(n, mode = "list")
+ for (i in 1:n) {
+ if (length(datalist[[i]]) == 1)
+ datalist.string[[i]] <- paste(names(datalist)[i],
+ "=", as.character(datalist[[i]]), sep = "")
+ if (is.vector(datalist[[i]]) && length(datalist[[i]]) > 1)
+ datalist.string[[i]] <- paste(names(datalist)[i],
+ "=c(", paste(as.character(datalist[[i]]), collapse = ", "),
+ ")", sep = "")
+ if (is.array(datalist[[i]]))
+ datalist.string[[i]] <- paste(names(datalist)[i],
+ "= structure(.Data= c(", paste(as.character(as.vector(aperm(datalist[[i]]))),
+ collapse = ", "), "), .Dim=c(", paste(as.character(dim(datalist[[i]])),
+ collapse = ", "), "))", sep = "")
+ }
+ datalist.tofile <- paste("list(", paste(unlist(datalist.string),
+ collapse = ", "), ")", sep = "")
+ return(datalist.tofile)
+}
Property changes on: trunk/R2WinBUGS/R/formatdata.R
___________________________________________________________________
Name: svn:keywords
+ Id
Added: trunk/R2WinBUGS/R/monitor.R
===================================================================
--- trunk/R2WinBUGS/R/monitor.R (rev 0)
+++ trunk/R2WinBUGS/R/monitor.R 2006-11-24 13:24:51 UTC (rev 3)
@@ -0,0 +1,51 @@
+"monitor" <-
+function (a, n.chains, trans=NULL, keep.all=FALSE, Rupper.keep=FALSE) {
+
+# If keep.all=T: a is a n x m x k array:
+# m sequences of length n, k variables measured
+# If keep.all=F: a is a 2n x m x k array (first half will be discarded)
+#
+# trans is a vector of length k: "" if no transformation, or "log" or "logit"
+# (If trans is not defined, it will be set to "log" for parameters that
+# are all-positive and 0 otherwise.)
+#
+# If Rupper.keep=T: keep Rupper. (Otherwise don't display it.)
+ invlogit <- function (x) {1 / (1 + exp(-x))}
+ nparams <- if(length(dim(a)) < 3) 1 else dim(a)[length(dim(a))]
+ # Calculation and initialization of the required matrix "output"
+ output <- matrix( , ncol = if(n.chains > 1){if(Rupper.keep) 10 else 9} else 7, nrow = nparams)
+ if (length(dim(a))==2) a <- array (a, c(dim(a),1))
+ if (!keep.all){
+ n <- floor(dim(a)[1]/2)
+ a <- a[(n+1):(2*n), , , drop = FALSE]
+ }
+ if (is.null(trans))
+ trans <- ifelse ((apply (a<=0, 3, sum))==0, "log", "")
+ for (i in 1:nparams){
+ # Rupper.keep: discard Rupper (nobody ever uses it)
+ ai <- a[ , , i, drop = FALSE]
+ if (trans[i]=="log"){
+ conv.p <- conv.par(log(ai), n.chains, Rupper.keep=Rupper.keep) # reason????
+ conv.p <- list(quantiles = exp(conv.p$quantiles),
+ confshrink = conv.p$confshrink, n.eff = conv.p$n.eff)
+ }
+ else if (trans[i]=="logit"){
+ conv.p <- conv.par(logit(ai), n.chains, Rupper.keep=Rupper.keep)
+ conv.p <- list(quantiles = invlogit(conv.p$quantiles),
+ confshrink = conv.p$confshrink, n.eff = conv.p$n.eff)
+ }
+ else conv.p <- conv.par(ai, n.chains, Rupper.keep=Rupper.keep)
+ output[i, ] <- c(mean(ai), sd(as.vector(ai)),
+ conv.p$quantiles,
+ if(n.chains > 1) conv.p$confshrink,
+ if(n.chains > 1) round(conv.p$n.eff, min(0, 1 - floor(log10(conv.p$n.eff))))
+ )
+ }
+ if(n.chains > 1)
+ dimnames(output) <- list(dimnames(a)[[3]], c("mean","sd",
+ "2.5%","25%","50%","75%","97.5%", "Rhat", if(Rupper.keep) "Rupper","n.eff"))
+ else
+ dimnames(output) <- list(dimnames(a)[[3]], c("mean","sd",
+ "2.5%","25%","50%","75%","97.5%"))
+ return (output)
+}
Property changes on: trunk/R2WinBUGS/R/monitor.R
___________________________________________________________________
Name: svn:keywords
+ Id
Added: trunk/R2WinBUGS/R/openbugs.R
===================================================================
--- trunk/R2WinBUGS/R/openbugs.R (rev 0)
+++ trunk/R2WinBUGS/R/openbugs.R 2006-11-24 13:24:51 UTC (rev 3)
@@ -0,0 +1,85 @@
+openbugs <- function(d...
[truncated message content] |
