R : Copyright 2005, The R Foundation for Statistical Computing Version 2.1.1 (2005-06-20), ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for a HTML browser interface to help. Type 'q()' to quit R. > ### *
> ### > attach(NULL, name = "CheckExEnv") > assign(".CheckExEnv", as.environment(2), pos = length(search())) # base > ## add some hooks to label plot pages for base and grid graphics > setHook("plot.new", ".newplot.hook") > setHook("persp", ".newplot.hook") > setHook("grid.newpage", ".gridplot.hook") > > assign("cleanEx", + function(env = .GlobalEnv) { + rm(list = ls(envir = env, all.names = TRUE), envir = env) + RNGkind("default", "default") + set.seed(1) + options(warn = 1) + delayedAssign("T", stop("T used instead of TRUE"), + assign.env = .CheckExEnv) + delayedAssign("F", stop("F used instead of FALSE"), + assign.env = .CheckExEnv) + sch <- search() + newitems <- sch[! sch %in% .oldSearch] + for(item in rev(newitems)) + eval(substitute(detach(item), list(item=item))) + missitems <- .oldSearch[! .oldSearch %in% sch] + if(length(missitems)) + warning("items ", paste(missitems, collapse=", "), + " have been removed from the search path") + }, + env = .CheckExEnv) > assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now > assign("ptime", proc.time(), env = .CheckExEnv) > grDevices::postscript("MCMCpack-Examples.ps") > assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv) > options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly")) > options(warn = 1) > library('MCMCpack') Loading required package: coda Loading required package: MASS ## ## Markov Chain Monte Carlo Package (MCMCpack) ## Copyright (C) 2003-2005 Andrew D. Martin and Kevin M. Quinn ## ## Support provided by the U.S. National Science Foundation ## (Grants SES-0350646 and SES-0350613) ## > > assign(".oldSearch", search(), env = .CheckExEnv) > assign(".oldNS", loadedNamespaces(), env = .CheckExEnv) > cleanEx(); ..nameEx <- "MCMCdynamicEI" > > ### * MCMCdynamicEI > > flush(stderr()); flush(stdout()) > > ### Name: MCMCdynamicEI > ### Title: Markov Chain Monte Carlo for Quinn's Dynamic Ecological > ### Inference Model > ### Aliases: MCMCdynamicEI > ### Keywords: models > > ### ** Examples > > ## Not run: > ##D ## simulated data example 1 > ##D set.seed(3920) > ##D n <- 100 > ##D r0 <- rpois(n, 2000) > ##D r1 <- round(runif(n, 100, 4000)) > ##D p0.true <- pnorm(-1.5 + 1:n/(n/2)) > ##D p1.true <- pnorm(1.0 - 1:n/(n/4)) > ##D y0 <- rbinom(n, r0, p0.true) > ##D y1 <- rbinom(n, r1, p1.true) > ##D c0 <- y0 + y1 > ##D c1 <- (r0+r1) - c0 > ##D > ##D ## plot data > ##D dtomogplot(r0, r1, c0, c1, delay=0.1) > ##D > ##D ## fit dynamic model > ##D post1 <- MCMCdynamicEI(r0,r1,c0,c1, mcmc=40000, thin=5, verbose=100, > ##D seed=list(NA, 1)) > ##D > ##D ## fit exchangeable hierarchical model > ##D post2 <- MCMChierEI(r0,r1,c0,c1, mcmc=40000, thin=5, verbose=100, > ##D seed=list(NA, 2)) > ##D > ##D p0meanDyn <- colMeans(post1)[1:n] > ##D p1meanDyn <- colMeans(post1)[(n+1):(2*n)] > ##D p0meanHier <- colMeans(post2)[1:n] > ##D p1meanHier <- colMeans(post2)[(n+1):(2*n)] > ##D > ##D ## plot truth and posterior means > ##D pairs(cbind(p0.true, p0meanDyn, p0meanHier, p1.true, p1meanDyn, p1meanHier)) > ##D > ##D ## simulated data example 2 > ##D set.seed(8722) > ##D n <- 100 > ##D r0 <- rpois(n, 2000) > ##D r1 <- round(runif(n, 100, 4000)) > ##D p0.true <- pnorm(-1.0 + sin(1:n/(n/4))) > ##D p1.true <- pnorm(0.0 - 2*cos(1:n/(n/9))) > ##D y0 <- rbinom(n, r0, p0.true) > ##D y1 <- rbinom(n, r1, p1.true) > ##D c0 <- y0 + y1 > ##D c1 <- (r0+r1) - c0 > ##D > ##D ## plot data > ##D dtomogplot(r0, r1, c0, c1, delay=0.1) > ##D > ##D ## fit dynamic model > ##D post1 <- MCMCdynamicEI(r0,r1,c0,c1, mcmc=40000, thin=5, verbose=100, > ##D seed=list(NA, 1)) > ##D > ##D ## fit exchangeable hierarchical model > ##D post2 <- MCMChierEI(r0,r1,c0,c1, mcmc=40000, thin=5, verbose=100, > ##D seed=list(NA, 2)) > ##D > ##D p0meanDyn <- colMeans(post1)[1:n] > ##D p1meanDyn <- colMeans(post1)[(n+1):(2*n)] > ##D p0meanHier <- colMeans(post2)[1:n] > ##D p1meanHier <- colMeans(post2)[(n+1):(2*n)] > ##D > ##D ## plot truth and posterior means > ##D pairs(cbind(p0.true, p0meanDyn, p0meanHier, p1.true, p1meanDyn, p1meanHier)) > ##D > ## End(Not run) > > > > cleanEx(); ..nameEx <- "MCMCfactanal" > > ### * MCMCfactanal > > flush(stderr()); flush(stdout()) > > ### Name: MCMCfactanal > ### Title: Markov Chain Monte Carlo for Normal Theory Factor Analysis Model > ### Aliases: MCMCfactanal > ### Keywords: models > > ### ** Examples > > ## Not run: > ##D ### An example using the formula interface > ##D data(swiss) > ##D posterior <- MCMCfactanal(~Agriculture+Examination+Education+Catholic > ##D +Infant.Mortality, factors=2, > ##D lambda.constraints=list(Examination=list(1,"+"), > ##D Examination=list(2,"-"), Education=c(2,0), > ##D Infant.Mortality=c(1,0)), > ##D verbose=0, store.scores=FALSE, a0=1, b0=0.15, > ##D data=swiss, burnin=5000, mcmc=50000, thin=20) > ##D plot(posterior) > ##D summary(posterior) > ##D > ##D ### An example using the matrix interface > ##D Y <- cbind(swiss$Agriculture, swiss$Examination, > ##D swiss$Education, swiss$Catholic, > ##D swiss$Infant.Mortality) > ##D colnames(Y) <- c("Agriculture", "Examination", "Education", "Catholic", > ##D "Infant.Mortality") > ##D post <- MCMCfactanal(Y, factors=2, > ##D lambda.constraints=list(Examination=list(1,"+"), > ##D Examination=list(2,"-"), Education=c(2,0), > ##D Infant.Mortality=c(1,0)), > ##D verbose=0, store.scores=FALSE, a0=1, b0=0.15, > ##D burnin=5000, mcmc=50000, thin=20) > ##D > ## End(Not run) > > > > cleanEx(); ..nameEx <- "MCMChierEI" > > ### * MCMChierEI > > flush(stderr()); flush(stdout()) > > ### Name: MCMChierEI > ### Title: Markov Chain Monte Carlo for Wakefield's Hierarchial Ecological > ### Inference Model > ### Aliases: MCMChierEI > ### Keywords: models > > ### ** Examples > > ## Not run: > ##D ## simulated data example > ##D set.seed(3920) > ##D n <- 100 > ##D r0 <- round(runif(n, 400, 1500)) > ##D r1 <- round(runif(n, 100, 4000)) > ##D p0.true <- pnorm(rnorm(n, m=0.5, s=0.25)) > ##D p1.true <- pnorm(rnorm(n, m=0.0, s=0.10)) > ##D y0 <- rbinom(n, r0, p0.true) > ##D y1 <- rbinom(n, r1, p1.true) > ##D c0 <- y0 + y1 > ##D c1 <- (r0+r1) - c0 > ##D > ##D ## plot data > ##D tomogplot(r0, r1, c0, c1) > ##D > ##D ## fit exchangeable hierarchical model > ##D post <- MCMChierEI(r0,r1,c0,c1, mcmc=40000, thin=5, verbose=100, > ##D seed=list(NA, 1)) > ##D > ##D p0meanHier <- colMeans(post)[1:n] > ##D p1meanHier <- colMeans(post)[(n+1):(2*n)] > ##D > ##D ## plot truth and posterior means > ##D pairs(cbind(p0.true, p0meanHier, p1.true, p1meanHier)) > ##D > ## End(Not run) > > > > cleanEx(); ..nameEx <- "MCMCirt1d" > > ### * MCMCirt1d > > flush(stderr()); flush(stdout()) > > ### Name: MCMCirt1d > ### Title: Markov Chain Monte Carlo for One Dimensional Item Response > ### Theory Model > ### Aliases: MCMCirt1d > ### Keywords: models > > ### ** Examples > > ## Not run: > ##D ## US Supreme Court Example with inequality constraints > ##D data(SupremeCourt) > ##D posterior1 <- MCMCirt1d(t(SupremeCourt), > ##D theta.constraints=list(Scalia="+", Ginsburg="-"), > ##D B0.alpha=.2, B0.beta=.2, > ##D burnin=500, mcmc=100000, thin=20, verbose=500, > ##D store.item=TRUE) > ##D geweke.diag(posterior1) > ##D plot(posterior1) > ##D summary(posterior1) > ##D > ##D ## US Senate Example with equality constraints > ##D data(Senate) > ##D Sen.rollcalls <- Senate[,6:677] > ##D posterior2 <- MCMCirt1d(Sen.rollcalls, > ##D theta.constraints=list(KENNEDY=-2, HELMS=2), > ##D burnin=2000, mcmc=100000, thin=20, verbose=500) > ##D geweke.diag(posterior2) > ##D plot(posterior2) > ##D summary(posterior2) > ##D > ## End(Not run) > > > > cleanEx(); ..nameEx <- "MCMCirtKd" > > ### * MCMCirtKd > > flush(stderr()); flush(stdout()) > > ### Name: MCMCirtKd > ### Title: Markov Chain Monte Carlo for K-Dimensional Item Response Theory > ### Model > ### Aliases: MCMCirtKd > ### Keywords: models > > ### ** Examples > > ## Not run: > ##D data(SupremeCourt) > ##D # note that the rownames (the item names) are "1", "2", etc > ##D posterior1 <- MCMCirtKd(t(SupremeCourt), dimensions=1, > ##D burnin=5000, mcmc=50000, thin=10, > ##D B0=.25, store.item=TRUE, > ##D item.constraints=list("1"=list(2,"-"))) > ##D plot(posterior1) > ##D summary(posterior1) > ##D > ##D data(Senate) > ##D Sen.rollcalls <- Senate[,6:677] > ##D posterior2 <- MCMCirtKd(Sen.rollcalls, dimensions=2, > ##D burnin=5000, mcmc=50000, thin=10, > ##D item.constraints=list(rc2=list(2,"-"), rc2=c(3,0), > ##D rc3=list(3,"-")), > ##D B0=.25) > ##D plot(posterior2) > ##D summary(posterior2) > ##D > ## End(Not run) > > > > cleanEx(); ..nameEx <- "MCMClogit" > > ### * MCMClogit > > flush(stderr()); flush(stdout()) > > ### Name: MCMClogit > ### Title: Markov Chain Monte Carlo for Logistic Regression > ### Aliases: MCMClogit > ### Keywords: models > > ### ** Examples > > ## Not run: > ##D data(birthwt) > ##D posterior <- MCMClogit(low~age+as.factor(race)+smoke, data=birthwt) > ##D plot(posterior) > ##D summary(posterior) > ##D > ## End(Not run) > > > > cleanEx(); ..nameEx <- "MCMCmetrop1R" > > ### * MCMCmetrop1R > > flush(stderr()); flush(stdout()) > > ### Name: MCMCmetrop1R > ### Title: Metropolis Sampling from User-Written R function > ### Aliases: MCMCmetrop1R > ### Keywords: models > > ### ** Examples > > ## Not run: > ##D > ##D ## logistic regression with an improper uniform prior > ##D ## X and y are passed as args to MCMCmetrop1R > ##D > ##D logitfun <- function(beta){ > ##D eta <- X %*% beta > ##D p <- 1.0/(1.0+exp(-eta)) > ##D sum( y * log(p) + (1-y)*log(1-p) ) > ##D } > ##D > ##D x1 <- rnorm(1000) > ##D x2 <- rnorm(1000) > ##D Xdata <- cbind(1,x1,x2) > ##D p <- exp(.5 - x1 + x2)/(1+exp(.5 - x1 + x2)) > ##D yvector <- rbinom(1000, 1, p) > ##D > ##D post.samp <- MCMCmetrop1R(logitfun, theta.init=c(0,0,0), > ##D X=Xdata, y=yvector, > ##D thin=1, mcmc=40000, burnin=500, > ##D tune=c(1.5, 1.5, 1.5), > ##D verbose=500, logfun=TRUE, optim.maxit=100) > ##D > ##D raftery.diag(post.samp) > ##D plot(post.samp) > ##D summary(post.samp) > ##D ## ################################################## > ##D > ##D > ##D ## logistic regression with an improper uniform prior > ##D ## X and y are now in the global environment > ##D > ##D logitfun <- function(beta){ > ##D eta <- X %*% beta > ##D p <- 1.0/(1.0+exp(-eta)) > ##D sum( y * log(p) + (1-y)*log(1-p) ) > ##D } > ##D > ##D x1 <- rnorm(1000) > ##D x2 <- rnorm(1000) > ##D X <- cbind(1,x1,x2) > ##D p <- exp(.5 - x1 + x2)/(1+exp(.5 - x1 + x2)) > ##D y <- rbinom(1000, 1, p) > ##D > ##D post.samp <- MCMCmetrop1R(logitfun, theta.init=c(0,0,0), > ##D thin=1, mcmc=40000, burnin=500, > ##D tune=c(1.5, 1.5, 1.5), > ##D verbose=500, logfun=TRUE, optim.maxit=100) > ##D > ##D raftery.diag(post.samp) > ##D plot(post.samp) > ##D summary(post.samp) > ##D ## ################################################## > ##D > ##D ## negative binomial regression with an improper unform prior > ##D ## X and y are passed as args to MCMCmetrop1R > ##D negbinfun <- function(theta){ > ##D k <- length(theta) > ##D beta <- theta[1:(k-1)] > ##D alpha <- exp(theta[k]) > ##D mu <- exp(X %*% beta) > ##D log.like <- sum( > ##D lgamma(y+alpha) - lfactorial(y) - lgamma(alpha) + > ##D alpha * log(alpha/(alpha+mu)) + > ##D y * log(mu/(alpha+mu)) > ##D ) > ##D } > ##D > ##D n <- 1000 > ##D x1 <- rnorm(n) > ##D x2 <- rnorm(n) > ##D XX <- cbind(1,x1,x2) > ##D mu <- exp(1.5+x1+2*x2)*rgamma(n,1) > ##D yy <- rpois(n, mu) > ##D > ##D post.samp <- MCMCmetrop1R(negbinfun, theta.init=c(0,0,0,0), y=yy, X=XX, > ##D thin=1, mcmc=35000, burnin=1000, > ##D tune=1.5, verbose=500, logfun=TRUE, > ##D optim.maxit=500, seed=list(NA,1)) > ##D raftery.diag(post.samp) > ##D plot(post.samp) > ##D summary(post.samp) > ##D ## ################################################## > ##D > ##D > ## End(Not run) > > > > cleanEx(); ..nameEx <- "MCMCmixfactanal" > > ### * MCMCmixfactanal > > flush(stderr()); flush(stdout()) > > ### Name: MCMCmixfactanal > ### Title: Markov Chain Monte Carlo for Mixed Data Factor Analysis Model > ### Aliases: MCMCmixfactanal > ### Keywords: models > > ### ** Examples > > ## Not run: > ##D data(PErisk) > ##D > ##D post <- MCMCmixfactanal(~courts+barb2+prsexp2+prscorr2+gdpw2, > ##D factors=1, data=PErisk, > ##D lambda.constraints = list(courts=list(2,"-")), > ##D burnin=5000, mcmc=1000000, thin=50, > ##D verbose=500, L0=.25, store.lambda=TRUE, > ##D store.scores=TRUE, tune=1.2) > ##D plot(post) > ##D summary(post) > ##D > ##D > ##D > ##D library(MASS) > ##D data(Cars93) > ##D attach(Cars93) > ##D new.cars <- data.frame(Price, MPG.city, MPG.highway, > ##D Cylinders, EngineSize, Horsepower, > ##D RPM, Length, Wheelbase, Width, Weight, Origin) > ##D rownames(new.cars) <- paste(Manufacturer, Model) > ##D detach(Cars93) > ##D > ##D # drop obs 57 (Mazda RX 7) b/c it has a rotary engine > ##D new.cars <- new.cars[-57,] > ##D # drop 3 cylinder cars > ##D new.cars <- new.cars[new.cars$Cylinders!=3,] > ##D # drop 5 cylinder cars > ##D new.cars <- new.cars[new.cars$Cylinders!=5,] > ##D > ##D new.cars$log.Price <- log(new.cars$Price) > ##D new.cars$log.MPG.city <- log(new.cars$MPG.city) > ##D new.cars$log.MPG.highway <- log(new.cars$MPG.highway) > ##D new.cars$log.EngineSize <- log(new.cars$EngineSize) > ##D new.cars$log.Horsepower <- log(new.cars$Horsepower) > ##D > ##D new.cars$Cylinders <- ordered(new.cars$Cylinders) > ##D new.cars$Origin <- ordered(new.cars$Origin) > ##D > ##D > ##D post <- MCMCmixfactanal(~log.Price+log.MPG.city+ > ##D log.MPG.highway+Cylinders+log.EngineSize+ > ##D log.Horsepower+RPM+Length+ > ##D Wheelbase+Width+Weight+Origin, data=new.cars, > ##D lambda.constraints=list(log.Horsepower=list(2,"+"), > ##D log.Horsepower=c(3,0), weight=list(3,"+")), > ##D factors=2, > ##D burnin=5000, mcmc=500000, thin=100, verbose=500, > ##D L0=.25, tune=3.0) > ##D plot(post) > ##D summary(post) > ##D > ## End(Not run) > > > > cleanEx(); ..nameEx <- "MCMCmnl" > > ### * MCMCmnl > > flush(stderr()); flush(stdout()) > > ### Name: MCMCmnl > ### Title: Markov Chain Monte Carlo for Multinomial Logistic Regression > ### Aliases: MCMCmnl > ### Keywords: models > > ### ** Examples > > ## Not run: > ##D data(Nethvote) > ##D > ##D ## just a choice-specific X var > ##D post1 <- MCMCmnl(vote ~ > ##D choicevar(distD66, "sqdist", "D66") + > ##D choicevar(distPvdA, "sqdist", "PvdA") + > ##D choicevar(distVVD, "sqdist", "VVD") + > ##D choicevar(distCDA, "sqdist", "CDA"), > ##D baseline="D66", mcmc.method="MH", B0=0, > ##D verbose=500, mcmc=100000, thin=10, tune=1.0, > ##D data=Nethvote) > ##D > ##D plot(post1) > ##D summary(post1) > ##D > ##D > ##D ## just individual-specific X vars > ##D post2<- MCMCmnl(vote ~ > ##D relig + class + income + educ + age + urban, > ##D baseline="D66", mcmc.method="MH", B0=0, > ##D verbose=500, mcmc=100000, thin=10, tune=0.5, > ##D data=Nethvote) > ##D > ##D plot(post2) > ##D summary(post2) > ##D > ##D > ##D ## both choice-specific and individual-specific X vars > ##D post3 <- MCMCmnl(vote ~ > ##D choicevar(distD66, "sqdist", "D66") + > ##D choicevar(distPvdA, "sqdist", "PvdA") + > ##D choicevar(distVVD, "sqdist", "VVD") + > ##D choicevar(distCDA, "sqdist", "CDA") + > ##D relig + class + income + educ + age + urban, > ##D baseline="D66", mcmc.method="MH", B0=0, > ##D verbose=500, mcmc=100000, thin=10, tune=0.5, > ##D data=Nethvote) > ##D > ##D plot(post3) > ##D summary(post3) > ##D > ##D ## numeric y variable > ##D nethvote$vote <- as.numeric(nethvote$vote) > ##D post4 <- MCMCmnl(vote ~ > ##D choicevar(distD66, "sqdist", "2") + > ##D choicevar(distPvdA, "sqdist", "3") + > ##D choicevar(distVVD, "sqdist", "4") + > ##D choicevar(distCDA, "sqdist", "1") + > ##D relig + class + income + educ + age + urban, > ##D baseline="2", mcmc.method="MH", B0=0, > ##D verbose=500, mcmc=100000, thin=10, tune=0.5, > ##D data=Nethvote) > ##D > ##D plot(post4) > ##D summary(post4) > ##D > ##D > ##D ## Simulated data example with nonconstant choiceset > ##D n <- 1000 > ##D y <- matrix(0, n, 4) > ##D colnames(y) <- c("a", "b", "c", "d") > ##D xa <- rnorm(n) > ##D xb <- rnorm(n) > ##D xc <- rnorm(n) > ##D xd <- rnorm(n) > ##D xchoice <- cbind(xa, xb, xc, xd) > ##D z <- rnorm(n) > ##D for (i in 1:n){ > ##D ## randomly determine choiceset (c is always in choiceset) > ##D choiceset <- c(3, sample(c(1,2,4), 2, replace=FALSE)) > ##D numer <- matrix(0, 4, 1) > ##D for (j in choiceset){ > ##D if (j == 3){ > ##D numer[j] <- exp(xchoice[i, j] ) > ##D } > ##D else { > ##D numer[j] <- exp(xchoice[i, j] - z[i] ) > ##D } > ##D } > ##D p <- numer / sum(numer) > ##D y[i,] <- rmultinom(1, 1, p) > ##D y[i,-choiceset] <- -999 > ##D } > ##D > ##D post5 <- MCMCmnl(y~choicevar(xa, "x", "a") + > ##D choicevar(xb, "x", "b") + > ##D choicevar(xc, "x", "c") + > ##D choicevar(xd, "x", "d") + z, > ##D baseline="c", verbose=500, > ##D mcmc=100000, thin=10, tune=.85) > ##D > ##D plot(post5) > ##D summary(post5) > ##D > ##D > ## End(Not run) > > > > cleanEx(); ..nameEx <- "MCMCoprobit" > > ### * MCMCoprobit > > flush(stderr()); flush(stdout()) > > ### Name: MCMCoprobit > ### Title: Markov Chain Monte Carlo for Ordered Probit Regression > ### Aliases: MCMCoprobit > ### Keywords: models > > ### ** Examples > > ## Not run: > ##D x1 <- rnorm(100); x2 <- rnorm(100); > ##D z <- 1.0 + x1*0.1 - x2*0.5 + rnorm(100); > ##D y <- z; y[z < 0] <- 0; y[z >= 0 & z < 1] <- 1; > ##D y[z >= 1 & z < 1.5] <- 2; y[z >= 1.5] <- 3; > ##D posterior <- MCMCoprobit(y ~ x1 + x2, tune=0.3, mcmc=20000) > ##D plot(posterior) > ##D summary(posterior) > ##D > ## End(Not run) > > > > cleanEx(); ..nameEx <- "MCMCordfactanal" > > ### * MCMCordfactanal > > flush(stderr()); flush(stdout()) > > ### Name: MCMCordfactanal > ### Title: Markov Chain Monte Carlo for Ordinal Data Factor Analysis Model > ### Aliases: MCMCordfactanal > ### Keywords: models > > ### ** Examples > > ## Not run: > ##D data(painters) > ##D new.painters <- painters[,1:4] > ##D cuts <- apply(new.painters, 2, quantile, c(.25, .50, .75)) > ##D for (i in 1:4){ > ##D new.painters[new.painters[,i] ##D new.painters[new.painters[,i] ##D new.painters[new.painters[,i] ##D new.painters[new.painters[,i]<100,i] <- 400 > ##D } > ##D > ##D posterior <- MCMCordfactanal(~Composition+Drawing+Colour+Expression, > ##D data=new.painters, factors=1, > ##D lambda.constraints=list(Drawing=list(2,"+")), > ##D burnin=5000, mcmc=500000, thin=200, verbose=500, > ##D L0=0.5, store.lambda=TRUE, > ##D store.scores=TRUE, tune=1.2) > ##D plot(posterior) > ##D summary(posterior) > ##D > ## End(Not run) > > > > cleanEx(); ..nameEx <- "MCMCpoisson" > > ### * MCMCpoisson > > flush(stderr()); flush(stdout()) > > ### Name: MCMCpoisson > ### Title: Markov Chain Monte Carlo for Poisson Regression > ### Aliases: MCMCpoisson > ### Keywords: models > > ### ** Examples > > ## Not run: > ##D counts <- c(18,17,15,20,10,20,25,13,12) > ##D outcome <- gl(3,1,9) > ##D treatment <- gl(3,3) > ##D posterior <- MCMCpoisson(counts ~ outcome + treatment) > ##D plot(posterior) > ##D summary(posterior) > ##D > ## End(Not run) > > > > cleanEx(); ..nameEx <- "MCMCprobit" > > ### * MCMCprobit > > flush(stderr()); flush(stdout()) > > ### Name: MCMCprobit > ### Title: Markov Chain Monte Carlo for Probit Regression > ### Aliases: MCMCprobit > ### Keywords: models > > ### ** Examples > > ## Not run: > ##D data(birthwt) > ##D posterior <- MCMCprobit(low~age+as.factor(race)+smoke, data=birthwt) > ##D plot(posterior) > ##D summary(posterior) > ##D > ## End(Not run) > > > > cleanEx(); ..nameEx <- "MCMCregress" > > ### * MCMCregress > > flush(stderr()); flush(stdout()) > > ### Name: MCMCregress > ### Title: Markov Chain Monte Carlo for Gaussian Linear Regression > ### Aliases: MCMCregress > ### Keywords: models > > ### ** Examples > > ## Not run: > ##D line <- list(X = c(-2,-1,0,1,2), Y = c(1,3,3,3,5)) > ##D posterior <- MCMCregress(Y~X, data=line, verbose=1000) > ##D plot(posterior) > ##D raftery.diag(posterior) > ##D summary(posterior) > ## End(Not run) > > > > cleanEx(); ..nameEx <- "MCMCtobit" > > ### * MCMCtobit > > flush(stderr()); flush(stdout()) > > ### Name: MCMCtobit > ### Title: Markov Chain Monte Carlo for Gaussian Linear Regression with a > ### Censored Dependent Variable > ### Aliases: MCMCtobit > ### Keywords: models > > ### ** Examples > > library(survival) Loading required package: splines > example(tobin) tobin> tfit <- survreg(Surv(durable, durable > 0, type = "left") ~ age + quant, data = tobin, dist = "gaussian") tobin> predict(tfit, type = "response") [1] -3.04968679 -4.31254182 -0.54163315 -0.25607164 -1.85017727 -2.40987803 [7] -3.50629220 -0.74041486 -4.05145594 -3.55880518 -0.32223237 -3.68044619 [13] -3.65997456 -2.63254564 0.22382063 0.02177674 -0.09571284 -3.17696755 [19] -0.61521215 -3.13913903 > summary(tfit) Call: survreg(formula = Surv(durable, durable > 0, type = "left") ~ age + quant, data = tobin, dist = "gaussian") Value Std. Error z p (Intercept) 15.1449 16.0795 0.942 3.46e-01 age -0.1291 0.2186 -0.590 5.55e-01 quant -0.0455 0.0583 -0.782 4.34e-01 Log(scale) 1.7179 0.3103 5.536 3.10e-08 Scale= 5.57 Gaussian distribution Loglik(model)= -28.9 Loglik(intercept only)= -29.5 Chisq= 1.1 on 2 degrees of freedom, p= 0.58 Number of Newton-Raphson Iterations: 3 n= 20 > tfit.mcmc <- MCMCtobit(durable ~ age + quant, data=tobin, mcmc=30000, + verbose=1000) MCMCtobit iteration 1 of 31000 beta = 14.53324 -0.17875 -0.02518 sigma2 = 8.41402 MCMCtobit iteration 1001 of 31000 beta = 52.55459 -0.55878 -0.15851 sigma2 = 454.21928 MCMCtobit iteration 2001 of 31000 beta = 19.96620 -0.17549 -0.07358 sigma2 = 73.64337 MCMCtobit iteration 3001 of 31000 beta = 215.85382 -4.63364 -0.39967 sigma2 = 7795.87892 MCMCtobit iteration 4001 of 31000 beta = 36.92052 -1.40173 0.05722 sigma2 = 818.75771 MCMCtobit iteration 5001 of 31000 beta = -20.32282 0.56543 -0.05620 sigma2 = 329.88211 MCMCtobit iteration 6001 of 31000 beta = -19.99355 0.29027 -0.02940 sigma2 = 158.42768 MCMCtobit iteration 7001 of 31000 beta = 54.44763 -0.19122 -0.19175 sigma2 = 84.79789 MCMCtobit iteration 8001 of 31000 beta = 21.27023 -1.47644 0.15638 sigma2 = 209.03197 MCMCtobit iteration 9001 of 31000 beta = 19.55942 0.00574 -0.11454 sigma2 = 178.02732 MCMCtobit iteration 10001 of 31000 beta = -23.13812 -0.37182 0.11873 sigma2 = 273.72350 MCMCtobit iteration 11001 of 31000 beta = -36.96670 -0.38160 0.19984 sigma2 = 230.54203 MCMCtobit iteration 12001 of 31000 beta = 27.92896 -0.05697 -0.12680 sigma2 = 54.93730 MCMCtobit iteration 13001 of 31000 beta = 38.62670 -0.31182 -0.11833 sigma2 = 84.74766 MCMCtobit iteration 14001 of 31000 beta = 20.04418 -0.19912 -0.05552 sigma2 = 32.96554 MCMCtobit iteration 15001 of 31000 beta = 4.07612 0.30036 -0.10512 sigma2 = 199.13058 MCMCtobit iteration 16001 of 31000 beta = 9.04425 0.12401 -0.08732 sigma2 = 49.02667 MCMCtobit iteration 17001 of 31000 beta = 23.73888 -0.07277 -0.09530 sigma2 = 61.10130 MCMCtobit iteration 18001 of 31000 beta = 11.13250 -0.18867 -0.02441 sigma2 = 21.97008 MCMCtobit iteration 19001 of 31000 beta = 50.45571 -0.33907 -0.21748 sigma2 = 282.25100 MCMCtobit iteration 20001 of 31000 beta = 43.70535 -0.17534 -0.18929 sigma2 = 91.96757 MCMCtobit iteration 21001 of 31000 beta = 0.78279 -0.01303 -0.00555 sigma2 = 38.34054 MCMCtobit iteration 22001 of 31000 beta = -8.33373 0.12648 0.00181 sigma2 = 38.81276 MCMCtobit iteration 23001 of 31000 beta = 23.37540 -0.54442 -0.02633 sigma2 = 84.00760 MCMCtobit iteration 24001 of 31000 beta = 2.21793 0.07881 -0.03422 sigma2 = 48.28074 MCMCtobit iteration 25001 of 31000 beta = 37.66846 0.03082 -0.20538 sigma2 = 136.72523 MCMCtobit iteration 26001 of 31000 beta = 36.23518 -0.06838 -0.14872 sigma2 = 97.39464 MCMCtobit iteration 27001 of 31000 beta = 23.71261 0.15860 -0.14811 sigma2 = 94.82360 MCMCtobit iteration 28001 of 31000 beta = 12.21582 0.10392 -0.09573 sigma2 = 271.14305 MCMCtobit iteration 29001 of 31000 beta = 27.27615 0.20916 -0.18059 sigma2 = 119.14210 MCMCtobit iteration 30001 of 31000 beta = 36.86630 -0.21344 -0.10971 sigma2 = 20.06024 > plot(tfit.mcmc) > raftery.diag(tfit.mcmc) Quantile (q) = 0.025 Accuracy (r) = +/- 0.005 Probability (s) = 0.95 Burn-in Total Lower bound Dependence (M) (N) (Nmin) factor (I) (Intercept) 12 17720 3746 4.73 age 24 26532 3746 7.08 quant 16 19660 3746 5.25 sigma2 15 16908 3746 4.51 > summary(tfit.mcmc) Iterations = 1001:31000 Thinning interval = 1 Number of chains = 1 Sample size per chain = 30000 1. Empirical mean and standard deviation for each variable, plus standard error of the mean: Mean SD Naive SE Time-series SE (Intercept) 18.54583 38.3732 0.2215477 0.324151 age -0.29027 0.5688 0.0032838 0.008482 quant -0.04656 0.1396 0.0008059 0.001004 sigma2 169.79694 357.0018 2.0611510 10.102583 2. Quantiles for each variable: 2.5% 25% 50% 75% 97.5% (Intercept) -54.1518 -0.3687 17.56604 36.21097 97.1546 age -1.5624 -0.5159 -0.22307 0.01334 0.6110 quant -0.3224 -0.1127 -0.04689 0.01877 0.2284 sigma2 19.7112 48.5186 86.89723 172.36508 790.9222 > > > > cleanEx(); ..nameEx <- "dirichlet" > > ### * dirichlet > > flush(stderr()); flush(stdout()) > > ### Name: Dirichlet > ### Title: The Dirichlet Distribution > ### Aliases: Dirichlet ddirichlet rdirichlet > ### Keywords: distribution > > ### ** Examples > > density <- ddirichlet(c(.1,.2,.7), c(1,1,1)) > draws <- rdirichlet(20, c(1,1,1) ) > > > > cleanEx(); ..nameEx <- "dtomog" > > ### * dtomog > > flush(stderr()); flush(stdout()) > > ### Name: dtomogplot > ### Title: Dynamic Tomography Plot > ### Aliases: dtomogplot > ### Keywords: hplot > > ### ** Examples > > ## Not run: > ##D ## simulated data example 1 > ##D set.seed(3920) > ##D n <- 100 > ##D r0 <- rpois(n, 2000) > ##D r1 <- round(runif(n, 100, 4000)) > ##D p0.true <- pnorm(-1.5 + 1:n/(n/2)) > ##D p1.true <- pnorm(1.0 - 1:n/(n/4)) > ##D y0 <- rbinom(n, r0, p0.true) > ##D y1 <- rbinom(n, r1, p1.true) > ##D c0 <- y0 + y1 > ##D c1 <- (r0+r1) - c0 > ##D > ##D ## plot data > ##D dtomogplot(r0, r1, c0, c1, delay=0.1) > ##D > ##D ## simulated data example 2 > ##D set.seed(8722) > ##D n <- 100 > ##D r0 <- rpois(n, 2000) > ##D r1 <- round(runif(n, 100, 4000)) > ##D p0.true <- pnorm(-1.0 + sin(1:n/(n/4))) > ##D p1.true <- pnorm(0.0 - 2*cos(1:n/(n/9))) > ##D y0 <- rbinom(n, r0, p0.true) > ##D y1 <- rbinom(n, r1, p1.true) > ##D c0 <- y0 + y1 > ##D c1 <- (r0+r1) - c0 > ##D > ##D ## plot data > ##D dtomogplot(r0, r1, c0, c1, delay=0.1) > ## End(Not run) > > > > cleanEx(); ..nameEx <- "invgamma" > > ### * invgamma > > flush(stderr()); flush(stdout()) > > ### Name: InvGamma > ### Title: The Inverse Gamma Distribution > ### Aliases: dinvgamma rinvgamma InvGamma > ### Keywords: distribution > > ### ** Examples > > density <- dinvgamma(4.2, 1.1) > draws <- rinvgamma(10, 3.2) > > > > cleanEx(); ..nameEx <- "iwishart" > > ### * iwishart > > flush(stderr()); flush(stdout()) > > ### Name: InvWishart > ### Title: The Inverse Wishart Distribution > ### Aliases: diwish riwish InvWishart > ### Keywords: distribution > > ### ** Examples > > density <- diwish(matrix(c(2,-.3,-.3,4),2,2), 3, matrix(c(1,.3,.3,1),2,2)) > draw <- riwish(3, matrix(c(1,.3,.3,1),2,2)) > > > > cleanEx(); ..nameEx <- "noncenhypergeom" > > ### * noncenhypergeom > > flush(stderr()); flush(stdout()) > > ### Name: NoncenHypergeom > ### Title: The Noncentral Hypergeometric Distribution > ### Aliases: NoncenHypergeom rnoncenhypergeom dnoncenhypergeom > ### Keywords: distribution > > ### ** Examples > > density <- dnoncenhypergeom(NA, 500, 500, 500, 6.0) > draws <- rnoncenhypergeom(10, 500, 500, 500, 6.0) > > > > cleanEx(); ..nameEx <- "readscythe" > > ### * readscythe > > flush(stderr()); flush(stdout()) > > ### Name: read.Scythe > ### Title: Read a Matrix from a File written by Scythe > ### Aliases: read.Scythe > ### Keywords: file > > ### ** Examples > > ## Not run: mymatrix <- read.Scythe("myfile.txt") > > > > cleanEx(); ..nameEx <- "tomog" > > ### * tomog > > flush(stderr()); flush(stdout()) > > ### Name: tomogplot > ### Title: Tomography Plot > ### Aliases: tomogplot > ### Keywords: hplot > > ### ** Examples > > r0 <- rpois(100, 500) > r1 <- rpois(100, 200) > c0 <- rpois(100, 100) > c1 <- (r0 + r1) - c0 > tomogplot(r0, r1, c0, c1) [1] 0 > > > > cleanEx(); ..nameEx <- "vech" > > ### * vech > > flush(stderr()); flush(stdout()) > > ### Name: vech > ### Title: Extract Lower Triangular Elements from a Symmetric Matrix > ### Aliases: vech > ### Keywords: manip > > ### ** Examples > > symmat <- matrix(c(1,2,3,4,2,4,5,6,3,5,7,8,4,6,8,9),4,4) > vech(symmat) [1] 1 2 3 4 4 5 6 7 8 9 > > > > cleanEx(); ..nameEx <- "wishart" > > ### * wishart > > flush(stderr()); flush(stdout()) > > ### Name: Wishart > ### Title: The Wishart Distribution > ### Aliases: dwish rwish Wishart > ### Keywords: distribution > > ### ** Examples > > density <- dwish(matrix(c(2,-.3,-.3,4),2,2), 3, matrix(c(1,.3,.3,1),2,2)) > draw <- rwish(3, matrix(c(1,.3,.3,1),2,2)) > > > > cleanEx(); ..nameEx <- "writescythe" > > ### * writescythe > > flush(stderr()); flush(stdout()) > > ### Name: write.Scythe > ### Title: Write a Matrix to a File to be Read by Scythe > ### Aliases: write.Scythe > ### Keywords: file > > ### ** Examples > > ## Not run: write.Scythe(mymatrix, "myfile.txt") > > > > cleanEx(); ..nameEx <- "xpnd" > > ### * xpnd > > flush(stderr()); flush(stdout()) > > ### Name: xpnd > ### Title: Expand a Vector into a Symmetric Matrix > ### Aliases: xpnd > ### Keywords: manip > > ### ** Examples > > xpnd(c(1,2,3,4,4,5,6,7,8,9),4) [,1] [,2] [,3] [,4] [1,] 1 2 3 4 [2,] 2 4 5 6 [3,] 3 5 7 8 [4,] 4 6 8 9 > > > > ### *