1
2 # test outfun (function)
3
4 epsilon <- 1e-15
5
6 library(mcmc)
7
8 RNGkind("Marsaglia-Multicarry")
9 set.seed(42)
10
11 n <- 100
12 rho <- 0.5
13 beta0 <- 0.25
14 beta1 <- 1
15 beta2 <- 0.5
16
17 x1 <- rnorm(n)
18 x2 <- rho * x1 + sqrt(1 - rho^2) * rnorm(n)
19 eta <- beta0 + beta1 * x1 + beta2 * x2
20 p <- 1 / (1 + exp(- eta))
21 y <- as.numeric(runif(n) < p)
22
23 out <- glm(y ~ x1 + x2, family = binomial())
24
25 logl <- function(beta) {
26     if (length(beta) != 3) stop("length(beta) != 3")
27     beta0 <- beta[1]
28     beta1 <- beta[2]
29     beta2 <- beta[3]
30     eta <- beta0 + beta1 * x1 + beta2 * x2
31     p <- exp(eta) / (1 + exp(eta))
32     return(sum(log(p[y == 1])) + sum(log(1 - p[y == 0])))
33 }
34
35 out.metro <- metrop(logl, coefficients(out), 1e3, scale = 0.01)
36 out.metro$accept
37
38 out.metro <- metrop(out.metro, scale = 0.1)
39 out.metro$accept
40
41 out.metro <- metrop(out.metro, scale = 0.5)
42 out.metro$accept
43
44 apply(out.metro$batch, 2, mean)
45
46 out.metro <- metrop(logl, as.numeric(coefficients(out)), 1e2,
47     scale = 0.5, debug = TRUE, outfun = function(x) c(x, x^2))
48
49 out.metro <- metrop(out.metro)
50 out.metro$outfun
51 dim(out.metro$batch)
52
53 logl <- function(beta, x1, x2, y) {
54     if (length(beta) != 3) stop("length(beta) != 3")
55     beta0 <- beta[1]
56     beta1 <- beta[2]
57     beta2 <- beta[3]
58     eta <- beta0 + beta1 * x1 + beta2 * x2
59     p <- exp(eta) / (1 + exp(eta))
60     return(sum(log(p[y == 1])) + sum(log(1 - p[y == 0])))
61 }
62
63 out.metro <- metrop(logl, as.numeric(coefficients(out)), 1e2,
64     scale = 0.5, debug = TRUE, x1 = x1, x2 = x2, y = y)
65 out.metro$lud
66 out.metro <- metrop(out.metro, x1 = x1, x2 = x2, y = y)
67 out.metro$lud
68
69