Nothing
### R code from vignette source 'morph.Rnw'
###################################################
### code chunk number 1: foo
###################################################
options(keep.source = TRUE, width = 60)
foo <- packageDescription("mcmc")
###################################################
### code chunk number 2: morph.Rnw:98-100
###################################################
library(mcmc)
h2 <- morph(b=1)
###################################################
### code chunk number 3: morph.Rnw:105-107
###################################################
lud <- function(x) dt(x, df=3, log=TRUE)
lud.induced <- h2$lud(lud)
###################################################
### code chunk number 4: morph.Rnw:110-114
###################################################
curve(exp(Vectorize(lud.induced)(x)), from = -3, to = 3, lty = 2,
xlab = "t", ylab = "density")
curve(exp(lud(x)), add = TRUE)
legend("topright", c("t density", "induced density"), lty=1:2)
###################################################
### code chunk number 5: morph.Rnw:123-129
###################################################
lud(1:4)
lud(1)
foo <- try(lud.induced(1:4))
class(foo)
cat(foo, "\n")
lud.induced(1)
###################################################
### code chunk number 6: set-seed
###################################################
set.seed(42)
###################################################
### code chunk number 7: morph.Rnw:146-147
###################################################
out <- morph.metrop(lud, 0, blen=100, nbatch=100, morph=morph(b=1))
###################################################
### code chunk number 8: morph.Rnw:153-155
###################################################
# adjust scale to find a roughly 20% acceptance rate
out$accept
###################################################
### code chunk number 9: morph.Rnw:161-163
###################################################
out <- morph.metrop(out, scale=4)
out$accept
###################################################
### code chunk number 10: fig0too
###################################################
acf(out$batch)
###################################################
### code chunk number 11: fig0
###################################################
acf(out$batch)
###################################################
### code chunk number 12: morph.Rnw:187-188
###################################################
t.test(out$batch)
###################################################
### code chunk number 13: morph.Rnw:191-193
###################################################
colMeans(out$batch)
apply(out$batch, 2, sd) / sqrt(out$nbatch)
###################################################
### code chunk number 14: unmorph-metrop-adjust
###################################################
out.unmorph <- metrop(lud, 0, blen=1000, nbatch=1)
out.unmorph$accept
out.unmorph <- metrop(out.unmorph, scale=4)
out.unmorph$accept
out.unmorph <- metrop(out.unmorph, scale=6)
out.unmorph$accept
###################################################
### code chunk number 15: unmorph-metrop-t-long-run
###################################################
lout <- suppressWarnings(try(load("morph1.rda"), silent = TRUE))
if (inherits(lout, "try-error")) {
out.unmorph <- metrop(out.unmorph, blen = 1e5, nbatch = 1e3)
save(out.unmorph, file = "morph1.rda")
} else {
.Random.seed <- out.unmorph$final.seed
}
out.unmorph$accept
###################################################
### code chunk number 16: fig4too
###################################################
foo <- as.vector(out.unmorph$batch)
qqnorm(foo)
qqline(foo)
###################################################
### code chunk number 17: fig4
###################################################
foo <- as.vector(out.unmorph$batch)
qqnorm(foo)
qqline(foo)
###################################################
### code chunk number 18: shapiro-wilk
###################################################
shapiro.test(foo)
###################################################
### code chunk number 19: morph-metrop-t-long-run
###################################################
lout <- suppressWarnings(try(load("morph2.rda"), silent = TRUE))
if (inherits(lout, "try-error")) {
out.morph <- morph.metrop(out, blen = 1e5, nbatch = 1e3)
save(out.morph, file = "morph2.rda")
} else {
.Random.seed <- out.morph$final.seed
}
out.morph$accept
###################################################
### code chunk number 20: fig5too
###################################################
foo <- as.vector(out.morph$batch)
qqnorm(foo)
qqline(foo)
###################################################
### code chunk number 21: fig5
###################################################
foo <- as.vector(out.morph$batch)
qqnorm(foo)
qqline(foo)
###################################################
### code chunk number 22: shapiro-wilk
###################################################
shapiro.test(foo)
###################################################
### code chunk number 23: def-posterior-binom
###################################################
lud.binom <- function(beta, M, x, n) {
MB <- M %*% beta
sum(x * MB) - sum(n * log(1 + exp(MB)))
}
###################################################
### code chunk number 24: convert
###################################################
dat <- as.data.frame(UCBAdmissions)
dat.split <- split(dat, dat$Admit)
dat.split <- lapply(dat.split,
function(d) {
val <- as.character(d$Admit[1])
d["Admit"] <- NULL
names(d)[names(d) == "Freq"] <- val
d
})
dat <- merge(dat.split[[1]], dat.split[[2]])
###################################################
### code chunk number 25: build-model-matrix
###################################################
formula <- cbind(Admitted, Rejected) ~ (Gender + Dept)^2
mf <- model.frame(formula, dat)
M <- model.matrix(formula, mf)
###################################################
### code chunk number 26: morph.Rnw:396-398
###################################################
xi <- 0.30
nu <- 5
###################################################
### code chunk number 27: lud-binom
###################################################
lud.berkeley <- function(B)
lud.binom(B, M, dat$Admitted + xi * nu, dat$Admitted + dat$Rejected + nu)
###################################################
### code chunk number 28: morph.Rnw:410-419
###################################################
berkeley.out <- morph.metrop(lud.berkeley, rep(0, ncol(M)), blen=1000,
nbatch=1, scale=0.1, morph=morph(p=3))
berkeley.out$accept
berkeley.out <- morph.metrop(berkeley.out, scale=0.05)
berkeley.out$accept
berkeley.out <- morph.metrop(berkeley.out, scale=0.02)
berkeley.out$accept
berkeley.out <- morph.metrop(berkeley.out, blen=10000)
berkeley.out$accept
###################################################
### code chunk number 29: morph.Rnw:422-423
###################################################
berkeley.out <- morph.metrop(berkeley.out, blen=1, nbatch=100000)
###################################################
### code chunk number 30: morph.Rnw:428-433
###################################################
beta <- setNames(colMeans(berkeley.out$batch), colnames(M))
MB <- M %*% beta
dat$p <- dat$Admitted / (dat$Admitted + dat$Rejected)
dat$p.post <- exp(MB) / (1 + exp(MB))
dat
###################################################
### code chunk number 31: calculate-posterior-probabilities
###################################################
posterior.probabilities <-
t(apply(berkeley.out$batch, 1,
function(r) {
eMB <- exp(M %*% r)
eMB / (1 + eMB)
}))
quants <- apply(posterior.probabilities, 2, quantile, prob=c(0.05, 0.95))
quants.str <- matrix(apply(quants, 2,
function(r) sprintf("[%0.2f, %0.2f]", r[1], r[2])),
nrow=2, byrow=TRUE)
###################################################
### code chunk number 32: fig1
###################################################
x <- (0:5) * 2 + 1
plot(x[c(1, 6)] + 0.5 * c(-1, 1), 0:1,
xlab="Department", ylab="Probability", xaxt="n", type="n")
axis(1, x, LETTERS[1:6])
for(i in 1:6) {
lines((x[i]-0.25)*c(1, 1), quants[1:2, i], lwd=2, col="gray")
lines((x[i] + 0.25) * c(1, 1), quants[1:2, i + 6], lwd=2, col="gray")
points(x[i] + 0.25 * c(-1, 1), dat$p.post[i + c(0, 6)], pch=c("F", "M"))
}
###################################################
### code chunk number 33: cauchy-data
###################################################
n <- 15
mu0 <- 50
sigma0 <- 10
x <- rcauchy(n, mu0, sigma0)
round(sort(x), 1)
###################################################
### code chunk number 34: cauchy-log-unnormalized-posterior
###################################################
lup <- function(theta) {
if (any(is.na(theta)))
stop("NA or NaN in input to log unnormalized density function")
mu <- theta[1]
sigma <- theta[2]
if (sigma <= 0) return(-Inf)
if (any(! is.finite(theta))) return(-Inf)
result <- sum(dcauchy(x, mu, sigma, log = TRUE)) - log(sigma)
if (! is.finite(result)) {
warning(paste("Oops! mu = ", mu, "and sigma =", sigma))
}
return(result)
}
###################################################
### code chunk number 35: cauchy-robust
###################################################
mu.twiddle <- median(x)
sigma.twiddle <- IQR(x)
c(mu.twiddle, sigma.twiddle)
###################################################
### code chunk number 36: cauchy-posterior-mode
###################################################
oout <- optim(c(mu.twiddle, sigma.twiddle), lup,
control = list(fnscale = -1), hessian = TRUE)
stopifnot(oout$convergence == 0)
mu.hat <- oout$par[1]
sigma.hat <- oout$par[2]
c(mu.hat, sigma.hat)
###################################################
### code chunk number 37: cauchy-hessian
###################################################
oout$hessian
###################################################
### code chunk number 38: cauchy-se
###################################################
sqrt(- 1 / diag(oout$hessian))
###################################################
### code chunk number 39: cauchy-doit
###################################################
moo <- morph(b = 0.5, r = 7, center = c(mu.hat, sigma.hat))
mout <- morph.metrop(lup, c(mu.hat, sigma.hat), 1e4,
scale = 3, morph = moo)
mout$accept
mout <- morph.metrop(mout)
###################################################
### code chunk number 40: cfig1too
###################################################
acf(mout$batch)
###################################################
### code chunk number 41: cfig1
###################################################
acf(mout$batch)
###################################################
### code chunk number 42: cfig2too
###################################################
mu <- mout$batch[ , 1]
i <- seq(1, mout$nbatch, by = 15)
out.sub <- density(mu[i])
out <- density(mu, bw = out.sub$bw)
plot(out)
###################################################
### code chunk number 43: cfig2
###################################################
mu <- mout$batch[ , 1]
i <- seq(1, mout$nbatch, by = 15)
out.sub <- density(mu[i])
out <- density(mu, bw = out.sub$bw)
plot(out)
###################################################
### code chunk number 44: cfig3
###################################################
sigma <- mout$batch[ , 2]
out.sub <- density(sigma[i])
out <- density(sigma, bw = out.sub$bw)
plot(out)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.