1 2R version 2.13.1 (2011-07-08) 3Copyright (C) 2011 The R Foundation for Statistical Computing 4ISBN 3-900051-07-0 5Platform: x86_64-pc-linux-gnu (64-bit) 6 7R is free software and comes with ABSOLUTELY NO WARRANTY. 8You are welcome to redistribute it under certain conditions. 9Type 'license()' or 'licence()' for distribution details. 10 11R is a collaborative project with many contributors. 12Type 'contributors()' for more information and 13'citation()' on how to cite R or R packages in publications. 14 15Type 'demo()' for some demos, 'help()' for on-line help, or 16'help.start()' for an HTML browser interface to help. 17Type 'q()' to quit R. 18 19> 20> # test outfun (function) 21> 22> epsilon <- 1e-15 23> 24> library(mcmc) 25> 26> RNGkind("Marsaglia-Multicarry") 27> set.seed(42) 28> 29> n <- 100 30> rho <- 0.5 31> beta0 <- 0.25 32> beta1 <- 1 33> beta2 <- 0.5 34> 35> x1 <- rnorm(n) 36> x2 <- rho * x1 + sqrt(1 - rho^2) * rnorm(n) 37> eta <- beta0 + beta1 * x1 + beta2 * x2 38> p <- 1 / (1 + exp(- eta)) 39> y <- as.numeric(runif(n) < p) 40> 41> out <- glm(y ~ x1 + x2, family = binomial()) 42> 43> logl <- function(beta) { 44+ if (length(beta) != 3) stop("length(beta) != 3") 45+ beta0 <- beta[1] 46+ beta1 <- beta[2] 47+ beta2 <- beta[3] 48+ eta <- beta0 + beta1 * x1 + beta2 * x2 49+ p <- exp(eta) / (1 + exp(eta)) 50+ return(sum(log(p[y == 1])) + sum(log(1 - p[y == 0]))) 51+ } 52> 53> out.metro <- metrop(logl, coefficients(out), 1e3, scale = 0.01) 54> out.metro$accept 55[1] 0.982 56> 57> out.metro <- metrop(out.metro, scale = 0.1) 58> out.metro$accept 59[1] 0.795 60> 61> out.metro <- metrop(out.metro, scale = 0.5) 62> out.metro$accept 63[1] 0.264 64> 65> apply(out.metro$batch, 2, mean) 66[1] 0.06080257 1.42304941 0.52634149 67> 68> out.metro <- metrop(logl, as.numeric(coefficients(out)), 1e2, 69+ scale = 0.5, debug = TRUE, outfun = function(x) c(x, x^2)) 70> 71> niter <- out.metro$nbatch * out.metro$blen * out.metro$nspac 72> niter == nrow(out.metro$current) 73[1] TRUE 74> niter == nrow(out.metro$proposal) 75[1] TRUE 76> all(out.metro$current[1, ] == out.metro$initial) 77[1] TRUE 78> all(out.metro$current[niter, ] == out.metro$final) | 79+ all(out.metro$proposal[niter, ] == out.metro$final) 80[1] TRUE 81> 82> .Random.seed <- out.metro$initial.seed 83> d <- ncol(out.metro$proposal) 84> n <- nrow(out.metro$proposal) 85> my.proposal <- matrix(NA, n, d) 86> my.u <- double(n) 87> ska <- out.metro$scale 88> for (i in 1:n) { 89+ my.proposal[i, ] <- out.metro$current[i, ] + ska * rnorm(d) 90+ if (is.na(out.metro$u[i])) { 91+ my.u[i] <- NA 92+ } else { 93+ my.u[i] <- runif(1) 94+ } 95+ } 96> max(abs(out.metro$proposal - my.proposal)) < epsilon 97[1] TRUE 98> all(is.na(out.metro$u) == is.na(my.u)) 99[1] TRUE 100> all(out.metro$u[!is.na(out.metro$u)] == my.u[!is.na(my.u)]) 101[1] TRUE 102> 103> my.curr.log.green <- apply(out.metro$current, 1, logl) 104> my.prop.log.green <- apply(out.metro$proposal, 1, logl) 105> all(is.na(out.metro$u) == (my.prop.log.green > my.curr.log.green)) 106[1] TRUE 107> foo <- my.prop.log.green - my.curr.log.green 108> max(abs(foo - out.metro$log.green)) < epsilon 109[1] TRUE 110> 111> my.accept <- is.na(my.u) | my.u < exp(foo) 112> sum(my.accept) == round(n * out.metro$accept) 113[1] TRUE 114> if (my.accept[niter]) { 115+ all(out.metro$proposal[niter, ] == out.metro$final) 116+ } else { 117+ all(out.metro$current[niter, ] == out.metro$final) 118+ } 119[1] TRUE 120> 121> my.current <- out.metro$current 122> my.current[my.accept, ] <- my.proposal[my.accept, ] 123> my.current <- rbind(out.metro$initial, my.current[- niter, ]) 124> max(abs(out.metro$current - my.current)) < epsilon 125[1] TRUE 126> 127> my.path <- matrix(NA, n, d) 128> my.path[my.accept, ] <- out.metro$proposal[my.accept, ] 129> my.path[! my.accept, ] <- out.metro$current[! my.accept, ] 130> nspac <- out.metro$nspac 131> 132> my.path <- my.path[seq(nspac, niter, by = nspac), ] 133> 134> fred <- t(apply(my.path, 1, out.metro$outfun)) 135> k <- ncol(fred) 136> 137> foom <- array(as.vector(t(fred)), c(k, out.metro$blen, out.metro$nbatch)) 138> boom <- t(apply(foom, c(1, 3), mean)) 139> 140> all(dim(boom) == dim(out.metro$batch)) 141[1] TRUE 142> max(abs(boom - out.metro$batch)) < epsilon 143[1] TRUE 144> 145> goom <- cbind(my.path, my.path^2) 146> all(dim(goom) == dim(out.metro$batch)) 147[1] TRUE 148> max(abs(goom - out.metro$batch)) < epsilon 149[1] TRUE 150> 151