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