R : Copyright 2003, The R Development Core Team
Version 1.6.2  (2003-01-10)

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.

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.

[Previously saved workspace restored]

> invisible(options(echo = TRUE))
> 
>  set.seed(42)
>  library(mcmc)
> 
>  foo <- read.table(url("http://www.stat.umn.edu/geyer/PhD/F03/logit.txt"),
+      header = TRUE)
> 
>  out <- glm(y ~ x1 + x2 + x3 + x4, data = foo, family = binomial())
> 
>  beta.init <- as.numeric(coefficients(out))
> 
>  x <- foo
>  x$y <- NULL
>  x <- as.matrix(x)
>  x <- cbind(1, x)
>  dimnames(x) <- NULL
> 
>  y <- foo$y
> 
>  upost <- function(beta, x, y) {
+      eta <- x %*% beta
+      p <- ifelse(eta > 0, 1 / (1 + exp(- eta)), exp(eta) / (1 + exp(eta)))
+      logl <- sum(log(p[y == 1])) + sum(log(1 - p[y == 0]))
+      return(logl + sum(dnorm(beta, 0, 2, log = TRUE)))
+  }
> 
>  out <- metropolis(upost, beta.init, 1e3, list(scaling = 0.1), x = x, y = y)
>  out$accept
[1] 0.755
> 
>  out <- metropolis(upost, out$final, 1e3, list(scaling = 0.2), x = x, y = y)
>  out$accept
[1] 0.541
> 
>  out <- metropolis(upost, out$final, 1e3, list(scaling = 0.3), x = x, y = y)
>  out$accept
[1] 0.347
> 
>  out <- metropolis(upost, out$final, 1e3, list(scaling = 0.4), x = x, y = y)
>  out$accept
[1] 0.239
> 
>  out <- metropolis(upost, out$final, 1e4, list(scaling = 0.4), x = x, y = y)
>  out$accept
[1] 0.2365
> 
>  # library(ts)
>  # acf(t(out$path))
>  #
>  # not shown but looks like lag 50 has negligible autocorrelation
> 
>  fred <- olbm(t(out$path), 50)
>  sqrt(diag(fred))
[1] 0.01263578 0.01387757 0.01502682 0.01449509 0.01594919
> 
>  out <- metropolis(upost, out$final, 3e4, list(scaling = 0.4), x = x, y = y)
>  out$accept
[1] 0.2346333
> 
>  apply(out$path, 1, mean)
[1] 0.6525116 0.8024205 1.1605730 0.5175203 0.7392830
>  fred <- olbm(t(out$path), 50)
>  sqrt(diag(fred))
[1] 0.006563992 0.008412268 0.008768438 0.008200941 0.009817598
> 
>  # o. k. for the means
>  # now for the variances
> 
>  mu <- apply(out$path, 1, mean)
>  sally <- sweep(out$path, 1, mu)
>  sigmasq <- apply(sally^2, 1, mean)
>  sigmasq 
[1] 0.09033471 0.12819881 0.13103139 0.12701575 0.15855383
>  fred <- olbm(t(sally^2), 50)
>  sqrt(diag(fred))
[1] 0.002383716 0.004035140 0.003881882 0.003457774 0.004525574
> 
>  #### also just go for it ####
> 
>  out <- metropolis(upost, out$final, 2e5, list(scaling = 0.4), x = x, y = y)
>  out$accept
[1] 0.23484
> 
>  mu <- apply(out$path, 1, mean)
>  mu
[1] 0.6583889 0.8008012 1.1705819 0.5016185 0.7268777
>  round(mu, 4)
[1] 0.6584 0.8008 1.1706 0.5016 0.7269
>  fred <- olbm(t(out$path), 50)
>  mu.mcse <- sqrt(diag(fred))
>  mu.mcse
[1] 0.002588477 0.003354069 0.003266676 0.003059423 0.003830972
>  round(mu.mcse, 4)
[1] 0.0026 0.0034 0.0033 0.0031 0.0038
> 
>  sally <- sweep(out$path, 1, mu)
>  sigmasq <- apply(sally^2, 1, mean)
>  sigmasq
[1] 0.09263352 0.13617168 0.13042531 0.12493805 0.16158943
>  round(sigmasq, 4)
[1] 0.0926 0.1362 0.1304 0.1249 0.1616
>  fred <- olbm(t(sally^2), 50)
>  sigmasq.mcse <- sqrt(diag(fred))
>  sigmasq.mcse
[1] 0.0009658537 0.0015694822 0.0014763757 0.0013095152 0.0018905816
>  round(sigmasq.mcse, 4)
[1] 0.0010 0.0016 0.0015 0.0013 0.0019
> 
>  sigma <- sqrt(sigmasq)
>  sigma.mcse <- sigmasq.mcse / (2 * sigma)
>  sigma
[1] 0.3043576 0.3690145 0.3611444 0.3534658 0.4019819
>  sigma.mcse
[1] 0.001586709 0.002126586 0.002044024 0.001852393 0.002351576
>  round(sigma, 4)
[1] 0.3044 0.3690 0.3611 0.3535 0.4020
>  round(sigma.mcse, 4)
[1] 0.0016 0.0021 0.0020 0.0019 0.0024
> 
>  system("cat /proc/cpuinfo")
processor	: 0
vendor_id	: AuthenticAMD
cpu family	: 6
model		: 10
model name	: AMD Athlon(tm) XP 3000+
stepping	: 0
cpu MHz		: 2162.740
cache size	: 512 KB
fdiv_bug	: no
hlt_bug		: no
f00f_bug	: no
coma_bug	: no
fpu		: yes
fpu_exception	: yes
cpuid level	: 1
wp		: yes
flags		: fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 mmx fxsr sse syscall mmxext 3dnowext 3dnow
bogomips	: 4312.26

> proc.time()
[1] 87.59  0.99 90.75  0.00  0.00
>