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