1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
|
### R code from vignette source 'MCMCintro.Rnw'
###################################################
### code chunk number 1: MCMCintro.Rnw:34-42
###################################################
minmaxpost <- function(theta, data){
mu <- theta[1]
sigma <- exp(theta[2])
dnorm(data$min, mu, sigma, log=TRUE) +
dnorm(data$max, mu, sigma, log=TRUE) +
(data$n - 2) * log(pnorm(data$max, mu, sigma) -
pnorm(data$min, mu, sigma))
}
###################################################
### code chunk number 2: MCMCintro.Rnw:51-55
###################################################
data <- list(n=10, min=52, max=84)
library(LearnBayes)
fit <- laplace(minmaxpost, c(70, 2), data)
fit
###################################################
### code chunk number 3: MCMCintro.Rnw:60-64
###################################################
mycontour(minmaxpost, c(45, 95, 1.5, 4), data,
xlab=expression(mu), ylab=expression(paste("log ",sigma)))
mycontour(lbinorm, c(45, 95, 1.5, 4),
list(m=fit$mode, v=fit$var), add=TRUE, col="red")
###################################################
### code chunk number 4: MCMCintro.Rnw:73-78
###################################################
mcmc.fit <- rwmetrop(minmaxpost,
list(var=fit$v, scale=3),
c(70, 2),
10000,
data)
###################################################
### code chunk number 5: MCMCintro.Rnw:82-83
###################################################
mcmc.fit$accept
###################################################
### code chunk number 6: MCMCintro.Rnw:88-92
###################################################
mycontour(minmaxpost, c(45, 95, 1.5, 4), data,
xlab=expression(mu),
ylab=expression(paste("log ",sigma)))
points(mcmc.fit$par)
###################################################
### code chunk number 7: MCMCintro.Rnw:105-110
###################################################
mu <- mcmc.fit$par[, 1]
sigma <- exp(mcmc.fit$par[, 2])
P.75 <- mu + 0.674 * sigma
plot(density(P.75),
main="Posterior Density of Upper Quartile")
|
