R/doingBayesianDataAnalysis.R
In ndbrown6/lazypuppy: Doing Bayesian Data Analysis 2nd Edition
####################
# CHAPTER 2
####################
'_Chapter2Figure2.5' <- function(y_prior_variance = 1000,
beta_prior_variance = 0.0001,
jags_n_chains = 3,
jags_n_adapt = 1000,
jags_n_iter = 1000,
coda_n_iter = 10000,
coda_n_thin = 10,
file_path = ".",
verbose = TRUE)
{
data('HtWtData110')
Y = HtWtData110 %>% .[["weight"]]
X = HtWtData110 %>% .[["height"]]
n = nrow(HtWtData110)
jags_ = paste0("model{
for (i in 1:n) {
Y[i] ~ dnorm(mu[i], sigma)
mu[i] <- beta[1] + beta[2]*X[i]
}
for (j in 1:2) {
beta[j] ~ dnorm(0, ", beta_prior_variance, ")
}
sigma ~ dunif(0, ", y_prior_variance, ")
}")
model_ = jags.model(file = textConnection(jags_),
data = list(Y = Y,
n = n,
X = X),
n.chains = jags_n_chains,
n.adapt = jags_n_adapt,
quiet = !verbose)
update(object = model_,
n.iter = jags_n_iter,
progress.bar = "text")
p_ = coda.samples(model = model_,
variable.names = c("beta", "sigma"),
n.iter = coda_n_iter,
thin = coda_n_thin,
progress.bar = "text")
beta_ = bind_cols(`beta_0` = as.vector(do.call(cbind, lapply(p_, function(x) { data.frame(x)$`beta.1.`}))),
`beta_1` = as.vector(do.call(cbind, lapply(p_, function(x) { data.frame(x)$`beta.2.`}))))
plot_ = HtWtData110 %>%
dplyr::mutate(facet = paste0("median slope = ", signif(median(beta_$`beta_1`),3),
" | median intercept = ", signif(median(beta_$`beta_0`),3))) %>%
ggplot(aes(x=height, y=weight)) +
geom_abline(data = beta_, aes(slope = `beta_1`, intercept = `beta_0`), color = "orange", size = .1) +
geom_point(shape = 19, size = 3, alpha = .5, color = "steelblue") +
theme_bw() +
xlab("\nHeight (inches)") +
ylab("Weight (lbs)\n") +
facet_wrap(~facet)
pdf(file = paste0(file_path, "/HtWtData110.pdf"), width=4, height=4)
print(plot_)
dev.off()
plot_ = beta_ %>%
dplyr::mutate(facet = paste0("median slope = ", signif(median(beta_$`beta_1`),3))) %>%
ggplot(aes(x=`beta_1`)) +
geom_histogram(bins = 30, fill = "salmon", color = "red", alpha = .5) +
theme_bw() +
xlab(expression(beta[1])) +
ylab ("Frequency\n") +
facet_wrap(~facet)
pdf(file = paste0(file_path, "/PosteriorBeta1.pdf"), width=4, height=4)
print(plot_)
dev.off()
return(invisible(p_))
}
####################
# CHAPTER 7
####################
'Chapter7Figure7.4' <- function(n = 100000,
mean = 0,
sd = .2,
N = 20,
z = 14,
a = 1,
b = 1,
file_path = NULL,
verbose = TRUE)
{
'likelihood' <- function(theta, N, z)
{
return(invisible((theta^z)*((1-theta)^(N-z))))
}
'prior' <- function(theta, a, b)
{
return(invisible(dbeta(x=theta, shape1=a, shape2=b)))
}
'proposal' <- function(theta, mean, sd)
{
return(invisible(max(theta+rnorm(1, mean=mean, sd=sd), 0)))
}
'pmove' <- function(theta_x, theta_y, N, z, a, b)
{
return(invisible(min(1, (likelihood(theta_y, N, z)*prior(theta_y, a, b))/(likelihood(theta_x, N, z)*prior(theta_x, a, b)))))
}
theta = vector(mode="numeric", length=n)
theta[1] = .1
if (verbose) {
pb = txtProgressBar(min=2, max=n, style=3)
}
for (i in 2:n) {
p_theta = proposal(theta[i-1], mean=mean, sd=sd)
p_move = pmove(theta_x = theta[i-1], theta_y = p_theta, N = N, z = z, a = a, b = b)
if (runif(1, 0, 1)<=p_move) {
theta[i] = p_theta
} else {
theta[i] = theta[i-1]
}
if (verbose) {
setTxtProgressBar(pb, i)
}
}
if (verbose) {
close(pb)
}
data_ = tibble(index = 1:n,
theta = theta)
hdi95 = HDIofMCMC(x=data_ %>% .[["theta"]], HDI=.95)
plot_x = data_ %>%
dplyr::mutate(facet = paste0("Posterior theta | proposal SD = ", sd)) %>%
ggplot(aes(x=`theta`)) +
geom_histogram(bins = 30, fill = "salmon", color = "red", alpha = .5) +
geom_segment(data = tibble(x = hdi95[1],
xend = hdi95[2],
y = 0,
yend = 0), aes(x=x, y=y, xend=xend, yend=yend), size=3) +
theme_bw() +
xlab(expression(theta)) +
ylab ("Frequency\n") +
xlim(c(0,1)) +
facet_wrap(~facet)
plot_y = data_ %>%
dplyr::filter(index <= 100) %>%
dplyr::mutate(facet = "Beginning of chain") %>%
ggplot(aes(x=`index`, y=`theta`)) +
geom_step(stat = "identity", direction = "hv", colour = "salmon", size = 1) +
facet_wrap(~facet) +
theme_bw() +
ylab(expression(theta)) +
xlab ("\nIteration\n") +
xlim(c(1,100)) +
ylim(c(0,1))
plot_z = data_ %>%
dplyr::filter(index >= (n-1000)) %>%
dplyr::mutate(facet = "End of chain") %>%
ggplot(aes(x=`index`, y=`theta`)) +
geom_step(stat = "identity", direction = "hv", colour = "salmon", size = 1) +
facet_wrap(~facet) +
theme_bw() +
ylab(expression(theta)) +
xlab ("\nIteration\n") +
xlim(c(n-1000,n)) +
ylim(c(0,1))
if (!is.null(file_path)) {
pdf(file = paste0(file_path, "/Posteriortheta.pdf"), width=4, height=4)
print(plot_x)
print(plot_y)
print(plot_z)
dev.off()
} else {
dev.new()
print(plot_x)
dev.new()
print(plot_y)
dev.new()
print(plot_z)
}
return(invisible(data_))
}
'Chapter7Figure7.5' <- function(n = 100,
a1 = 2,
b1 = 2,
a2 = 2,
b2 = 2,
z1 = 6,
N1 = 8,
z2 = 2,
N2 = 7,
file_path = NULL)
{
'prior' <- function(theta1, theta2, a1, b1, a2, b2)
{
return(invisible(dbeta(theta1, a1, b1)*dbeta(theta2, a2, b2)))
}
'likelihood' <- function(theta1, theta2, N1, z1, N2, z2)
{
return(invisible((theta1^z1)*(1-theta1)^(N1-z1) * (theta2^z2)*(1-theta2)^(N2-z2)))
}
'posterior' <- function(theta1, theta2, a1, a2, b1, b2, N1, z1, N2, z2)
{
return(invisible(dbeta(theta1, a1+z1, b1+N1-z1)*dbeta(theta2, a2+z2, b2+N2-z2)))
}
theta1 = theta2 = seq(from=0, to=1, length=n)
ptheta = matrix(NA, nrow=n, ncol=n)
for (i in 1:n) {
for (j in 1:n) {
ptheta[i,j] = prior(theta1[i], theta2[j], a1, b1, a2, b2)
}
}
ll = matrix(NA, nrow=n, ncol=n)
for (i in 1:n) {
for (j in 1:n) {
ll[i,j] = likelihood(theta1[i], theta2[j], N1, z1, N2, z2)
}
}
pthetaD = matrix(NA, nrow=n, ncol=n)
for (i in 1:n) {
for (j in 1:n) {
pthetaD[i,j] = posterior(theta1[i], theta2[j], a1, a2, b1, b2, N1, z1, N2, z2)
}
}
if (!is.null(file_path)) {
pdf(file = paste0(file_path, "/Posteriortheta.pdf"), width=7, height=11)
par(mfrow=c(3,2))
persp(ptheta, xlim=c(0,1), ylim=c(0,1), zlim=c(0,5),
xlab = "\ntheta[1]",
ylab = "\ntheta[2]",
zlab = "\np(theta[1],theta[2])",
main = "Prior",
theta=35, phi=20, axes=TRUE, col="steelblue", shade=0.001, border="skyblue", ticktype="detailed")
contour(ptheta, xlim=c(0,1), ylim=c(0,1),
xlab = expression(theta[1]),
ylab = expression(theta[2]),
main = expression("p("~theta[1]~","~theta[2]~")"),
col="steelblue", las=1)
persp(ll, xlim=c(0,1), ylim=c(0,1),
xlab = "\ntheta[1]",
ylab = "\ntheta[2]",
zlab = "\np(D|theta[1],theta[2])",
main = "Likelihood",
theta=35, phi=20, axes=TRUE, col="steelblue", shade=0.001, border="skyblue", ticktype="detailed")
contour(ll, xlim=c(0,1), ylim=c(0,1),
xlab = expression(theta[1]),
ylab = expression(theta[2]),
main = expression("p(D|"~theta[1]~","~theta[2]~")"),
col="steelblue", las=1)
persp(pthetaD, xlim=c(0,1), ylim=c(0,1),
xlab = "\ntheta[1]",
ylab = "\ntheta[2]",
zlab = "\np(theta[1],theta[2]|D)",
main = "Posterior",
theta=35, phi=20, axes=TRUE, col="steelblue", shade=0.001, border="skyblue", ticktype="detailed")
contour(pthetaD, xlim=c(0,1), ylim=c(0,1),
xlab = expression(theta[1]),
ylab = expression(theta[2]),
main = expression("p("~theta[1]~","~theta[2]~"|D)"),
col="steelblue", las=1)
dev.off()
} else {
dev.new()
par(mfrow=c(1,2))
persp(ptheta, xlim=c(0,1), ylim=c(0,1), zlim=c(0,5),
xlab = "\ntheta[1]",
ylab = "\ntheta[2]",
zlab = "\np(theta[1],theta[2])",
main = "Prior",
theta=35, phi=20, axes=TRUE, col="steelblue", shade=0.001, border="skyblue", ticktype="detailed")
contour(ptheta, xlim=c(0,1), ylim=c(0,1),
xlab = expression(theta[1]),
ylab = expression(theta[2]),
main = expression("p("~theta[1]~","~theta[2]~")"),
col="steelblue", las=1)
dev.new()
par(mfrow=c(1,2))
persp(ll, xlim=c(0,1), ylim=c(0,1),
xlab = "\ntheta[1]",
ylab = "\ntheta[2]",
zlab = "\np(D|theta[1],theta[2])",
main = "Likelihood",
theta=35, phi=20, axes=TRUE, col="steelblue", shade=0.001, border="skyblue", ticktype="detailed")
contour(ll, xlim=c(0,1), ylim=c(0,1),
xlab = expression(theta[1]),
ylab = expression(theta[2]),
main = expression("p(D|"~theta[1]~","~theta[2]~")"),
col="steelblue", las=1)
dev.new()
par(mfrow=c(1,2))
persp(pthetaD, xlim=c(0,1), ylim=c(0,1),
xlab = "\ntheta[1]",
ylab = "\ntheta[2]",
zlab = "\np(theta[1],theta[2]|D)",
main = "Posterior",
theta=35, phi=20, axes=TRUE, col="steelblue", shade=0.001, border="skyblue", ticktype="detailed")
contour(pthetaD, xlim=c(0,1), ylim=c(0,1),
xlab = expression(theta[1]),
ylab = expression(theta[2]),
main = expression("p("~theta[1]~","~theta[2]~"|D)"),
col="steelblue", las=1)
}
return(invisible(list(prior = ptheta,
likelihood = ll,
posterior = pthetaD)))
}
####################
# CHAPTER 9
####################
'Chapter9Figure9.2' <- function(theta = seq(from=0, to=1, length=100),
omega = seq(from=0, to=1, length=100),
K = 100,
A_omega = 2,
B_omega = 2,
N = 12,
z = 9,
file_path = NULL)
{
'pthetaomega' <- function(theta, omega, K, A_omega, B_omega)
{
p_omega = dbeta(x=omega, shape1=A_omega, shape2=B_omega)
p_thetaomega = dbeta(x=theta, shape1=(omega*(K-2)+1), shape2=((1-omega)*(K-2))+1)
return(invisible(p_omega*p_thetaomega))
}
prior_matrix = matrix(NA, length(theta), length(omega))
for (i in 1:length(theta)) {
for (j in 1:length(omega)) {
prior_matrix[i,j] = pthetaomega(theta=theta[i], omega=omega[j], K=K, A_omega=A_omega, B_omega=B_omega)
}
}
prior_matrix = (prior_matrix/sum(prior_matrix))/(0.01^2)
'pomega' <- function(omega, A_omega, B_omega)
{
p_omega = dbeta(omega, shape1=A_omega, shape2=B_omega)
return(invisible(p_omega))
}
p_omega = pomega(omega=omega, A_omega=A_omega, B_omega=B_omega)
'pthetaKnomega' <- function(theta, omega, K, A_omega, B_omega)
{
p_theta = pthetaomega(theta, omega, K, A_omega, B_omega)
return(invisible(p_theta))
}
p_theta = apply(prior_matrix, 1, sum)
p_theta = (p_theta/sum(p_theta))/(0.01)
p_thetao_.25 = pthetaKnomega(theta, omega=.25, K=K, A_omega=A_omega, B_omega=B_omega)
p_thetao_.25 = (p_thetao_.25/sum(p_thetao_.25))/(0.01)
p_thetao_.75 = pthetaKnomega(theta, omega=.75, K=K, A_omega=A_omega, B_omega=B_omega)
p_thetao_.75 = (p_thetao_.75/sum(p_thetao_.75))/(0.01)
if (!is.null(file_path)) {
pdf(file=paste0(file_path, "/Prior.pdf"), width=8, height=6)
par(mfrow=c(2,3))
persp(prior_matrix, xlim=c(0,1), ylim=c(0,1),
xlab = "\ntheta",
ylab = "\nomega",
zlab = "\n\n\np(theta,omega)",
main = "Prior",
theta=-25, phi=30, axes=TRUE, col="steelblue", shade=0.001, border="skyblue", ticktype="detailed")
contour(prior_matrix, xlim=c(0,1), ylim=c(0,1),
xlab = expression(theta),
ylab = expression(omega),
main = expression("p("~theta~"|"~omega~")"),
col="steelblue", las=1)
plot(omega, p_omega, type="l", col="steelblue",
xlab = expression(omega),
ylab = expression(p(omega)),
main = expression("Marginal"~p(omega)),
las = 1, lwd=2)
plot(theta, p_theta, type="l", col="steelblue",
xlab = expression(theta),
ylab = expression(p(theta)),
main = expression("Marginal"~p(theta)),
las = 1, lwd=2)
plot(theta, p_thetao_.25, type="l", col="steelblue",
xlab = expression(theta),
ylab = expression(p(theta~"|"~omega~"=.25")),
main = expression("Marginal"~p(theta~"|"~omega~"=.25")),
las = 1, lwd=2)
plot(theta, p_thetao_.75, type="l", col="steelblue",
xlab = expression(theta),
ylab = expression(p(theta~"|"~omega~"=.75")),
main = expression("Marginal"~p(theta~"|"~omega~"=.75")),
las = 1, lwd=2)
dev.off()
} else {
dev.new()
par(mfrow=c(2,3))
persp(prior_matrix, xlim=c(0,1), ylim=c(0,1),
xlab = "\ntheta",
ylab = "\nomega",
zlab = "\n\n\np(theta,omega)",
main = "Prior",
theta=-25, phi=30, axes=TRUE, col="steelblue", shade=0.001, border="skyblue", ticktype="detailed")
contour(prior_matrix, xlim=c(0,1), ylim=c(0,1),
xlab = expression(theta),
ylab = expression(omega),
main = expression("p("~theta~"|"~omega~")"),
col="steelblue", las=1)
plot(omega, p_omega, type="l", col="steelblue",
xlab = expression(omega),
ylab = expression(p(omega)),
main = expression("Marginal"~p(omega)),
las = 1, lwd=2)
plot(theta, p_theta, type="l", col="steelblue",
xlab = expression(theta),
ylab = expression(p(theta)),
main = expression("Marginal"~p(theta)),
las = 1, lwd=2)
plot(theta, p_thetao_.25, type="l", col="steelblue",
xlab = expression(theta),
ylab = expression(p(theta~"|"~omega~"=.25")),
main = expression("Marginal"~p(theta~"|"~omega~"=.25")),
las = 1, lwd=2)
plot(theta, p_thetao_.75, type="l", col="steelblue",
xlab = expression(theta),
ylab = expression(p(theta~"|"~omega~"=.75")),
main = expression("Marginal"~p(theta~"|"~omega~"=.75")),
las = 1, lwd=2)
}
'likelihood' <- function(theta, N, z)
{
return(invisible((theta^z)*((1-theta)^(N-z))))
}
likelihood_matrix = matrix(NA, length(theta), length(omega))
for (i in 1:length(theta)) {
for (j in 1:length(omega)) {
likelihood_matrix[i,j] = likelihood(theta=theta[i], N=N, z=z)
}
}
if (!is.null(file_path)) {
pdf(file=paste0(file_path, "/Likelihood.pdf"), width=8, height=5)
par(mfrow=c(1,2))
persp(likelihood_matrix, xlim=c(0,1), ylim=c(0,1),
xlab = "\ntheta",
ylab = "\nomega",
zlab = "\n\n\np(D|theta,omega)",
main = "Likelihood",
theta=-25, phi=30, axes=TRUE, col="steelblue", shade=0.001, border="skyblue", ticktype="detailed")
contour(likelihood_matrix, xlim=c(0,1), ylim=c(0,1),
xlab = expression(theta),
ylab = expression(omega),
main = expression("p("~D~"|"~theta~","~omega~")"),
col="steelblue", las=1)
dev.off()
} else {
dev.new()
par(mfrow=c(1,2))
persp(likelihood_matrix, xlim=c(0,1), ylim=c(0,1),
xlab = "\ntheta",
ylab = "\nomega",
zlab = "\n\n\np(D|theta,omega)",
main = "Likelihood",
theta=-25, phi=30, axes=TRUE, col="steelblue", shade=0.001, border="skyblue", ticktype="detailed")
contour(likelihood_matrix, xlim=c(0,1), ylim=c(0,1),
xlab = expression(theta),
ylab = expression(omega),
main = expression("p("~D~"|"~theta~","~omega~")"),
col="steelblue", las=1)
}
'posterior' <- function(theta, omega, K, A_omega, B_omega, N, z)
{
prior_matrix = matrix(NA, length(theta), length(omega))
for (i in 1:length(theta)) {
for (j in 1:length(omega)) {
prior_matrix[i,j] = pthetaomega(theta=theta[i], omega=omega[j], K=K, A_omega=A_omega, B_omega=B_omega)
}
}
prior_matrix = (prior_matrix/sum(prior_matrix))/(0.01^2)
likelihood_matrix = matrix(NA, length(theta), length(omega))
for (i in 1:length(theta)) {
for (j in 1:length(omega)) {
likelihood_matrix[i,j] = likelihood(theta=theta[i], N=N, z=z)
}
}
posterior_matrix = prior_matrix*likelihood_matrix
posterior_matrix = posterior_matrix/sum(posterior_matrix)
return(invisible(posterior_matrix))
}
posterior_matrix = posterior(theta=theta, omega=omega, K=K, A_omega=A_omega, B_omega=B_omega, N=N, z=z)
posterior_matrix = posterior_matrix/(0.01^2)
p_omegaD = apply(posterior_matrix, 2, sum)
p_omegaD = (p_omegaD/sum(p_omegaD))/(0.01)
p_thetaD = apply(posterior_matrix, 1, sum)
p_thetaD = (p_thetaD/sum(p_thetaD))/(0.01)
p_thetaD_.25 = posterior(theta, omega=.25, K=K, A_omega=A_omega, B_omega=B_omega, N, z)
p_thetaD_.25 = (p_thetaD_.25/sum(p_thetaD_.25))/(0.01)
p_thetaD_.75 = posterior(theta, omega=.75, K=K, A_omega=A_omega, B_omega=B_omega, N, z)
p_thetaD_.75 = (p_thetaD_.75/sum(p_thetaD_.75))/(0.01)
if (!is.null(file_path)) {
pdf(file=paste0(file_path, "/Posterior.pdf"), width=8, height=6)
par(mfrow=c(2,3))
persp(posterior_matrix, xlim=c(0,1), ylim=c(0,1),
xlab = "\ntheta",
ylab = "\nomega",
zlab = "\n\n\np(theta,omega|D)",
main = "Posterior",
theta=-25, phi=30, axes=TRUE, col="steelblue", shade=0.001, border="skyblue", ticktype="detailed")
contour(posterior_matrix, xlim=c(0,1), ylim=c(0,1),
xlab = expression(theta),
ylab = expression(omega),
main = expression("p("~theta~","~omega~"|D)"),
col="steelblue", las=1)
plot(omega, p_omegaD, type="l", col="steelblue",
xlab = expression(omega),
ylab = expression("p("~omega~"|D)"),
main = expression("Marginal p("~omega~"|D)"),
las = 1, lwd=2)
plot(omega, p_thetaD, type="l", col="steelblue",
xlab = expression(theta),
ylab = expression("p("~theta~"|D)"),
main = expression("Marginal p("~theta~"|D)"),
las = 1, lwd=2)
plot(theta, p_thetaD_.25, type="l", col="steelblue",
xlab = expression(theta),
ylab = expression(p(theta~","~omega~"=.25|D")),
main = expression("Marginal"~p(theta~","~omega~"=.25|D")),
las = 1, lwd=2)
plot(theta, p_thetaD_.75, type="l", col="steelblue",
xlab = expression(theta),
ylab = expression(p(theta~","~omega~"=.75|D")),
main = expression("Marginal"~p(theta~","~omega~"=.75|D")),
las = 1, lwd=2)
dev.off()
} else {
dev.new()
par(mfrow=c(2,3))
persp(posterior_matrix, xlim=c(0,1), ylim=c(0,1),
xlab = "\ntheta",
ylab = "\nomega",
zlab = "\n\n\np(theta,omega|D)",
main = "Posterior",
theta=-25, phi=30, axes=TRUE, col="steelblue", shade=0.001, border="skyblue", ticktype="detailed")
contour(posterior_matrix, xlim=c(0,1), ylim=c(0,1),
xlab = expression(theta),
ylab = expression(omega),
main = expression("p("~theta~","~omega~"|D)"),
col="steelblue", las=1)
plot(omega, p_omegaD, type="l", col="steelblue",
xlab = expression(omega),
ylab = expression("p("~omega~"|D)"),
main = expression("Marginal p("~omega~"|D)"),
las = 1, lwd=2)
plot(omega, p_thetaD, type="l", col="steelblue",
xlab = expression(theta),
ylab = expression("p("~theta~"|D)"),
main = expression("Marginal p("~theta~"|D)"),
las = 1, lwd=2)
plot(theta, p_thetaD_.25, type="l", col="steelblue",
xlab = expression(theta),
ylab = expression(p(theta~","~omega~"=.25|D")),
main = expression("Marginal"~p(theta~","~omega~"=.25|D")),
las = 1, lwd=2)
plot(theta, p_thetaD_.75, type="l", col="steelblue",
xlab = expression(theta),
ylab = expression(p(theta~","~omega~"=.75|D")),
main = expression("Marginal"~p(theta~","~omega~"=.75|D")),
las = 1, lwd=2)
}
return(invisible(list(prior = prior_matrix,
likelihood = likelihood_matrix,
posterior = posterior_matrix)))
}
'Chapter9Figure9.3' <- function(theta = seq(from=0, to=1, length=100),
omega = seq(from=0, to=1, length=100),
K = 6,
A_omega = 20,
B_omega = 20,
N = 12,
z = 9,
file_path = NULL)
{
p_ = Chapter9Figure9.2(K = K,
A_omega = A_omega,
B_omega = B_omega,
N = N,
z = z,
file_path = file_path)
return(invisible(p_))
}
ndbrown6/lazypuppy documentation built on Sept. 14, 2020, 7:33 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
####################
# CHAPTER 2
####################
'_Chapter2Figure2.5' <- function(y_prior_variance = 1000,
beta_prior_variance = 0.0001,
jags_n_chains = 3,
jags_n_adapt = 1000,
jags_n_iter = 1000,
coda_n_iter = 10000,
coda_n_thin = 10,
file_path = ".",
verbose = TRUE)
{
data('HtWtData110')
Y = HtWtData110 %>% .[["weight"]]
X = HtWtData110 %>% .[["height"]]
n = nrow(HtWtData110)
jags_ = paste0("model{
for (i in 1:n) {
Y[i] ~ dnorm(mu[i], sigma)
mu[i] <- beta[1] + beta[2]*X[i]
}
for (j in 1:2) {
beta[j] ~ dnorm(0, ", beta_prior_variance, ")
}
sigma ~ dunif(0, ", y_prior_variance, ")
}")
model_ = jags.model(file = textConnection(jags_),
data = list(Y = Y,
n = n,
X = X),
n.chains = jags_n_chains,
n.adapt = jags_n_adapt,
quiet = !verbose)
update(object = model_,
n.iter = jags_n_iter,
progress.bar = "text")
p_ = coda.samples(model = model_,
variable.names = c("beta", "sigma"),
n.iter = coda_n_iter,
thin = coda_n_thin,
progress.bar = "text")
beta_ = bind_cols(`beta_0` = as.vector(do.call(cbind, lapply(p_, function(x) { data.frame(x)$`beta.1.`}))),
`beta_1` = as.vector(do.call(cbind, lapply(p_, function(x) { data.frame(x)$`beta.2.`}))))
plot_ = HtWtData110 %>%
dplyr::mutate(facet = paste0("median slope = ", signif(median(beta_$`beta_1`),3),
" | median intercept = ", signif(median(beta_$`beta_0`),3))) %>%
ggplot(aes(x=height, y=weight)) +
geom_abline(data = beta_, aes(slope = `beta_1`, intercept = `beta_0`), color = "orange", size = .1) +
geom_point(shape = 19, size = 3, alpha = .5, color = "steelblue") +
theme_bw() +
xlab("\nHeight (inches)") +
ylab("Weight (lbs)\n") +
facet_wrap(~facet)
pdf(file = paste0(file_path, "/HtWtData110.pdf"), width=4, height=4)
print(plot_)
dev.off()
plot_ = beta_ %>%
dplyr::mutate(facet = paste0("median slope = ", signif(median(beta_$`beta_1`),3))) %>%
ggplot(aes(x=`beta_1`)) +
geom_histogram(bins = 30, fill = "salmon", color = "red", alpha = .5) +
theme_bw() +
xlab(expression(beta[1])) +
ylab ("Frequency\n") +
facet_wrap(~facet)
pdf(file = paste0(file_path, "/PosteriorBeta1.pdf"), width=4, height=4)
print(plot_)
dev.off()
return(invisible(p_))
}
####################
# CHAPTER 7
####################
'Chapter7Figure7.4' <- function(n = 100000,
mean = 0,
sd = .2,
N = 20,
z = 14,
a = 1,
b = 1,
file_path = NULL,
verbose = TRUE)
{
'likelihood' <- function(theta, N, z)
{
return(invisible((theta^z)*((1-theta)^(N-z))))
}
'prior' <- function(theta, a, b)
{
return(invisible(dbeta(x=theta, shape1=a, shape2=b)))
}
'proposal' <- function(theta, mean, sd)
{
return(invisible(max(theta+rnorm(1, mean=mean, sd=sd), 0)))
}
'pmove' <- function(theta_x, theta_y, N, z, a, b)
{
return(invisible(min(1, (likelihood(theta_y, N, z)*prior(theta_y, a, b))/(likelihood(theta_x, N, z)*prior(theta_x, a, b)))))
}
theta = vector(mode="numeric", length=n)
theta[1] = .1
if (verbose) {
pb = txtProgressBar(min=2, max=n, style=3)
}
for (i in 2:n) {
p_theta = proposal(theta[i-1], mean=mean, sd=sd)
p_move = pmove(theta_x = theta[i-1], theta_y = p_theta, N = N, z = z, a = a, b = b)
if (runif(1, 0, 1)<=p_move) {
theta[i] = p_theta
} else {
theta[i] = theta[i-1]
}
if (verbose) {
setTxtProgressBar(pb, i)
}
}
if (verbose) {
close(pb)
}
data_ = tibble(index = 1:n,
theta = theta)
hdi95 = HDIofMCMC(x=data_ %>% .[["theta"]], HDI=.95)
plot_x = data_ %>%
dplyr::mutate(facet = paste0("Posterior theta | proposal SD = ", sd)) %>%
ggplot(aes(x=`theta`)) +
geom_histogram(bins = 30, fill = "salmon", color = "red", alpha = .5) +
geom_segment(data = tibble(x = hdi95[1],
xend = hdi95[2],
y = 0,
yend = 0), aes(x=x, y=y, xend=xend, yend=yend), size=3) +
theme_bw() +
xlab(expression(theta)) +
ylab ("Frequency\n") +
xlim(c(0,1)) +
facet_wrap(~facet)
plot_y = data_ %>%
dplyr::filter(index <= 100) %>%
dplyr::mutate(facet = "Beginning of chain") %>%
ggplot(aes(x=`index`, y=`theta`)) +
geom_step(stat = "identity", direction = "hv", colour = "salmon", size = 1) +
facet_wrap(~facet) +
theme_bw() +
ylab(expression(theta)) +
xlab ("\nIteration\n") +
xlim(c(1,100)) +
ylim(c(0,1))
plot_z = data_ %>%
dplyr::filter(index >= (n-1000)) %>%
dplyr::mutate(facet = "End of chain") %>%
ggplot(aes(x=`index`, y=`theta`)) +
geom_step(stat = "identity", direction = "hv", colour = "salmon", size = 1) +
facet_wrap(~facet) +
theme_bw() +
ylab(expression(theta)) +
xlab ("\nIteration\n") +
xlim(c(n-1000,n)) +
ylim(c(0,1))
if (!is.null(file_path)) {
pdf(file = paste0(file_path, "/Posteriortheta.pdf"), width=4, height=4)
print(plot_x)
print(plot_y)
print(plot_z)
dev.off()
} else {
dev.new()
print(plot_x)
dev.new()
print(plot_y)
dev.new()
print(plot_z)
}
return(invisible(data_))
}
'Chapter7Figure7.5' <- function(n = 100,
a1 = 2,
b1 = 2,
a2 = 2,
b2 = 2,
z1 = 6,
N1 = 8,
z2 = 2,
N2 = 7,
file_path = NULL)
{
'prior' <- function(theta1, theta2, a1, b1, a2, b2)
{
return(invisible(dbeta(theta1, a1, b1)*dbeta(theta2, a2, b2)))
}
'likelihood' <- function(theta1, theta2, N1, z1, N2, z2)
{
return(invisible((theta1^z1)*(1-theta1)^(N1-z1) * (theta2^z2)*(1-theta2)^(N2-z2)))
}
'posterior' <- function(theta1, theta2, a1, a2, b1, b2, N1, z1, N2, z2)
{
return(invisible(dbeta(theta1, a1+z1, b1+N1-z1)*dbeta(theta2, a2+z2, b2+N2-z2)))
}
theta1 = theta2 = seq(from=0, to=1, length=n)
ptheta = matrix(NA, nrow=n, ncol=n)
for (i in 1:n) {
for (j in 1:n) {
ptheta[i,j] = prior(theta1[i], theta2[j], a1, b1, a2, b2)
}
}
ll = matrix(NA, nrow=n, ncol=n)
for (i in 1:n) {
for (j in 1:n) {
ll[i,j] = likelihood(theta1[i], theta2[j], N1, z1, N2, z2)
}
}
pthetaD = matrix(NA, nrow=n, ncol=n)
for (i in 1:n) {
for (j in 1:n) {
pthetaD[i,j] = posterior(theta1[i], theta2[j], a1, a2, b1, b2, N1, z1, N2, z2)
}
}
if (!is.null(file_path)) {
pdf(file = paste0(file_path, "/Posteriortheta.pdf"), width=7, height=11)
par(mfrow=c(3,2))
persp(ptheta, xlim=c(0,1), ylim=c(0,1), zlim=c(0,5),
xlab = "\ntheta[1]",
ylab = "\ntheta[2]",
zlab = "\np(theta[1],theta[2])",
main = "Prior",
theta=35, phi=20, axes=TRUE, col="steelblue", shade=0.001, border="skyblue", ticktype="detailed")
contour(ptheta, xlim=c(0,1), ylim=c(0,1),
xlab = expression(theta[1]),
ylab = expression(theta[2]),
main = expression("p("~theta[1]~","~theta[2]~")"),
col="steelblue", las=1)
persp(ll, xlim=c(0,1), ylim=c(0,1),
xlab = "\ntheta[1]",
ylab = "\ntheta[2]",
zlab = "\np(D|theta[1],theta[2])",
main = "Likelihood",
theta=35, phi=20, axes=TRUE, col="steelblue", shade=0.001, border="skyblue", ticktype="detailed")
contour(ll, xlim=c(0,1), ylim=c(0,1),
xlab = expression(theta[1]),
ylab = expression(theta[2]),
main = expression("p(D|"~theta[1]~","~theta[2]~")"),
col="steelblue", las=1)
persp(pthetaD, xlim=c(0,1), ylim=c(0,1),
xlab = "\ntheta[1]",
ylab = "\ntheta[2]",
zlab = "\np(theta[1],theta[2]|D)",
main = "Posterior",
theta=35, phi=20, axes=TRUE, col="steelblue", shade=0.001, border="skyblue", ticktype="detailed")
contour(pthetaD, xlim=c(0,1), ylim=c(0,1),
xlab = expression(theta[1]),
ylab = expression(theta[2]),
main = expression("p("~theta[1]~","~theta[2]~"|D)"),
col="steelblue", las=1)
dev.off()
} else {
dev.new()
par(mfrow=c(1,2))
persp(ptheta, xlim=c(0,1), ylim=c(0,1), zlim=c(0,5),
xlab = "\ntheta[1]",
ylab = "\ntheta[2]",
zlab = "\np(theta[1],theta[2])",
main = "Prior",
theta=35, phi=20, axes=TRUE, col="steelblue", shade=0.001, border="skyblue", ticktype="detailed")
contour(ptheta, xlim=c(0,1), ylim=c(0,1),
xlab = expression(theta[1]),
ylab = expression(theta[2]),
main = expression("p("~theta[1]~","~theta[2]~")"),
col="steelblue", las=1)
dev.new()
par(mfrow=c(1,2))
persp(ll, xlim=c(0,1), ylim=c(0,1),
xlab = "\ntheta[1]",
ylab = "\ntheta[2]",
zlab = "\np(D|theta[1],theta[2])",
main = "Likelihood",
theta=35, phi=20, axes=TRUE, col="steelblue", shade=0.001, border="skyblue", ticktype="detailed")
contour(ll, xlim=c(0,1), ylim=c(0,1),
xlab = expression(theta[1]),
ylab = expression(theta[2]),
main = expression("p(D|"~theta[1]~","~theta[2]~")"),
col="steelblue", las=1)
dev.new()
par(mfrow=c(1,2))
persp(pthetaD, xlim=c(0,1), ylim=c(0,1),
xlab = "\ntheta[1]",
ylab = "\ntheta[2]",
zlab = "\np(theta[1],theta[2]|D)",
main = "Posterior",
theta=35, phi=20, axes=TRUE, col="steelblue", shade=0.001, border="skyblue", ticktype="detailed")
contour(pthetaD, xlim=c(0,1), ylim=c(0,1),
xlab = expression(theta[1]),
ylab = expression(theta[2]),
main = expression("p("~theta[1]~","~theta[2]~"|D)"),
col="steelblue", las=1)
}
return(invisible(list(prior = ptheta,
likelihood = ll,
posterior = pthetaD)))
}
####################
# CHAPTER 9
####################
'Chapter9Figure9.2' <- function(theta = seq(from=0, to=1, length=100),
omega = seq(from=0, to=1, length=100),
K = 100,
A_omega = 2,
B_omega = 2,
N = 12,
z = 9,
file_path = NULL)
{
'pthetaomega' <- function(theta, omega, K, A_omega, B_omega)
{
p_omega = dbeta(x=omega, shape1=A_omega, shape2=B_omega)
p_thetaomega = dbeta(x=theta, shape1=(omega*(K-2)+1), shape2=((1-omega)*(K-2))+1)
return(invisible(p_omega*p_thetaomega))
}
prior_matrix = matrix(NA, length(theta), length(omega))
for (i in 1:length(theta)) {
for (j in 1:length(omega)) {
prior_matrix[i,j] = pthetaomega(theta=theta[i], omega=omega[j], K=K, A_omega=A_omega, B_omega=B_omega)
}
}
prior_matrix = (prior_matrix/sum(prior_matrix))/(0.01^2)
'pomega' <- function(omega, A_omega, B_omega)
{
p_omega = dbeta(omega, shape1=A_omega, shape2=B_omega)
return(invisible(p_omega))
}
p_omega = pomega(omega=omega, A_omega=A_omega, B_omega=B_omega)
'pthetaKnomega' <- function(theta, omega, K, A_omega, B_omega)
{
p_theta = pthetaomega(theta, omega, K, A_omega, B_omega)
return(invisible(p_theta))
}
p_theta = apply(prior_matrix, 1, sum)
p_theta = (p_theta/sum(p_theta))/(0.01)
p_thetao_.25 = pthetaKnomega(theta, omega=.25, K=K, A_omega=A_omega, B_omega=B_omega)
p_thetao_.25 = (p_thetao_.25/sum(p_thetao_.25))/(0.01)
p_thetao_.75 = pthetaKnomega(theta, omega=.75, K=K, A_omega=A_omega, B_omega=B_omega)
p_thetao_.75 = (p_thetao_.75/sum(p_thetao_.75))/(0.01)
if (!is.null(file_path)) {
pdf(file=paste0(file_path, "/Prior.pdf"), width=8, height=6)
par(mfrow=c(2,3))
persp(prior_matrix, xlim=c(0,1), ylim=c(0,1),
xlab = "\ntheta",
ylab = "\nomega",
zlab = "\n\n\np(theta,omega)",
main = "Prior",
theta=-25, phi=30, axes=TRUE, col="steelblue", shade=0.001, border="skyblue", ticktype="detailed")
contour(prior_matrix, xlim=c(0,1), ylim=c(0,1),
xlab = expression(theta),
ylab = expression(omega),
main = expression("p("~theta~"|"~omega~")"),
col="steelblue", las=1)
plot(omega, p_omega, type="l", col="steelblue",
xlab = expression(omega),
ylab = expression(p(omega)),
main = expression("Marginal"~p(omega)),
las = 1, lwd=2)
plot(theta, p_theta, type="l", col="steelblue",
xlab = expression(theta),
ylab = expression(p(theta)),
main = expression("Marginal"~p(theta)),
las = 1, lwd=2)
plot(theta, p_thetao_.25, type="l", col="steelblue",
xlab = expression(theta),
ylab = expression(p(theta~"|"~omega~"=.25")),
main = expression("Marginal"~p(theta~"|"~omega~"=.25")),
las = 1, lwd=2)
plot(theta, p_thetao_.75, type="l", col="steelblue",
xlab = expression(theta),
ylab = expression(p(theta~"|"~omega~"=.75")),
main = expression("Marginal"~p(theta~"|"~omega~"=.75")),
las = 1, lwd=2)
dev.off()
} else {
dev.new()
par(mfrow=c(2,3))
persp(prior_matrix, xlim=c(0,1), ylim=c(0,1),
xlab = "\ntheta",
ylab = "\nomega",
zlab = "\n\n\np(theta,omega)",
main = "Prior",
theta=-25, phi=30, axes=TRUE, col="steelblue", shade=0.001, border="skyblue", ticktype="detailed")
contour(prior_matrix, xlim=c(0,1), ylim=c(0,1),
xlab = expression(theta),
ylab = expression(omega),
main = expression("p("~theta~"|"~omega~")"),
col="steelblue", las=1)
plot(omega, p_omega, type="l", col="steelblue",
xlab = expression(omega),
ylab = expression(p(omega)),
main = expression("Marginal"~p(omega)),
las = 1, lwd=2)
plot(theta, p_theta, type="l", col="steelblue",
xlab = expression(theta),
ylab = expression(p(theta)),
main = expression("Marginal"~p(theta)),
las = 1, lwd=2)
plot(theta, p_thetao_.25, type="l", col="steelblue",
xlab = expression(theta),
ylab = expression(p(theta~"|"~omega~"=.25")),
main = expression("Marginal"~p(theta~"|"~omega~"=.25")),
las = 1, lwd=2)
plot(theta, p_thetao_.75, type="l", col="steelblue",
xlab = expression(theta),
ylab = expression(p(theta~"|"~omega~"=.75")),
main = expression("Marginal"~p(theta~"|"~omega~"=.75")),
las = 1, lwd=2)
}
'likelihood' <- function(theta, N, z)
{
return(invisible((theta^z)*((1-theta)^(N-z))))
}
likelihood_matrix = matrix(NA, length(theta), length(omega))
for (i in 1:length(theta)) {
for (j in 1:length(omega)) {
likelihood_matrix[i,j] = likelihood(theta=theta[i], N=N, z=z)
}
}
if (!is.null(file_path)) {
pdf(file=paste0(file_path, "/Likelihood.pdf"), width=8, height=5)
par(mfrow=c(1,2))
persp(likelihood_matrix, xlim=c(0,1), ylim=c(0,1),
xlab = "\ntheta",
ylab = "\nomega",
zlab = "\n\n\np(D|theta,omega)",
main = "Likelihood",
theta=-25, phi=30, axes=TRUE, col="steelblue", shade=0.001, border="skyblue", ticktype="detailed")
contour(likelihood_matrix, xlim=c(0,1), ylim=c(0,1),
xlab = expression(theta),
ylab = expression(omega),
main = expression("p("~D~"|"~theta~","~omega~")"),
col="steelblue", las=1)
dev.off()
} else {
dev.new()
par(mfrow=c(1,2))
persp(likelihood_matrix, xlim=c(0,1), ylim=c(0,1),
xlab = "\ntheta",
ylab = "\nomega",
zlab = "\n\n\np(D|theta,omega)",
main = "Likelihood",
theta=-25, phi=30, axes=TRUE, col="steelblue", shade=0.001, border="skyblue", ticktype="detailed")
contour(likelihood_matrix, xlim=c(0,1), ylim=c(0,1),
xlab = expression(theta),
ylab = expression(omega),
main = expression("p("~D~"|"~theta~","~omega~")"),
col="steelblue", las=1)
}
'posterior' <- function(theta, omega, K, A_omega, B_omega, N, z)
{
prior_matrix = matrix(NA, length(theta), length(omega))
for (i in 1:length(theta)) {
for (j in 1:length(omega)) {
prior_matrix[i,j] = pthetaomega(theta=theta[i], omega=omega[j], K=K, A_omega=A_omega, B_omega=B_omega)
}
}
prior_matrix = (prior_matrix/sum(prior_matrix))/(0.01^2)
likelihood_matrix = matrix(NA, length(theta), length(omega))
for (i in 1:length(theta)) {
for (j in 1:length(omega)) {
likelihood_matrix[i,j] = likelihood(theta=theta[i], N=N, z=z)
}
}
posterior_matrix = prior_matrix*likelihood_matrix
posterior_matrix = posterior_matrix/sum(posterior_matrix)
return(invisible(posterior_matrix))
}
posterior_matrix = posterior(theta=theta, omega=omega, K=K, A_omega=A_omega, B_omega=B_omega, N=N, z=z)
posterior_matrix = posterior_matrix/(0.01^2)
p_omegaD = apply(posterior_matrix, 2, sum)
p_omegaD = (p_omegaD/sum(p_omegaD))/(0.01)
p_thetaD = apply(posterior_matrix, 1, sum)
p_thetaD = (p_thetaD/sum(p_thetaD))/(0.01)
p_thetaD_.25 = posterior(theta, omega=.25, K=K, A_omega=A_omega, B_omega=B_omega, N, z)
p_thetaD_.25 = (p_thetaD_.25/sum(p_thetaD_.25))/(0.01)
p_thetaD_.75 = posterior(theta, omega=.75, K=K, A_omega=A_omega, B_omega=B_omega, N, z)
p_thetaD_.75 = (p_thetaD_.75/sum(p_thetaD_.75))/(0.01)
if (!is.null(file_path)) {
pdf(file=paste0(file_path, "/Posterior.pdf"), width=8, height=6)
par(mfrow=c(2,3))
persp(posterior_matrix, xlim=c(0,1), ylim=c(0,1),
xlab = "\ntheta",
ylab = "\nomega",
zlab = "\n\n\np(theta,omega|D)",
main = "Posterior",
theta=-25, phi=30, axes=TRUE, col="steelblue", shade=0.001, border="skyblue", ticktype="detailed")
contour(posterior_matrix, xlim=c(0,1), ylim=c(0,1),
xlab = expression(theta),
ylab = expression(omega),
main = expression("p("~theta~","~omega~"|D)"),
col="steelblue", las=1)
plot(omega, p_omegaD, type="l", col="steelblue",
xlab = expression(omega),
ylab = expression("p("~omega~"|D)"),
main = expression("Marginal p("~omega~"|D)"),
las = 1, lwd=2)
plot(omega, p_thetaD, type="l", col="steelblue",
xlab = expression(theta),
ylab = expression("p("~theta~"|D)"),
main = expression("Marginal p("~theta~"|D)"),
las = 1, lwd=2)
plot(theta, p_thetaD_.25, type="l", col="steelblue",
xlab = expression(theta),
ylab = expression(p(theta~","~omega~"=.25|D")),
main = expression("Marginal"~p(theta~","~omega~"=.25|D")),
las = 1, lwd=2)
plot(theta, p_thetaD_.75, type="l", col="steelblue",
xlab = expression(theta),
ylab = expression(p(theta~","~omega~"=.75|D")),
main = expression("Marginal"~p(theta~","~omega~"=.75|D")),
las = 1, lwd=2)
dev.off()
} else {
dev.new()
par(mfrow=c(2,3))
persp(posterior_matrix, xlim=c(0,1), ylim=c(0,1),
xlab = "\ntheta",
ylab = "\nomega",
zlab = "\n\n\np(theta,omega|D)",
main = "Posterior",
theta=-25, phi=30, axes=TRUE, col="steelblue", shade=0.001, border="skyblue", ticktype="detailed")
contour(posterior_matrix, xlim=c(0,1), ylim=c(0,1),
xlab = expression(theta),
ylab = expression(omega),
main = expression("p("~theta~","~omega~"|D)"),
col="steelblue", las=1)
plot(omega, p_omegaD, type="l", col="steelblue",
xlab = expression(omega),
ylab = expression("p("~omega~"|D)"),
main = expression("Marginal p("~omega~"|D)"),
las = 1, lwd=2)
plot(omega, p_thetaD, type="l", col="steelblue",
xlab = expression(theta),
ylab = expression("p("~theta~"|D)"),
main = expression("Marginal p("~theta~"|D)"),
las = 1, lwd=2)
plot(theta, p_thetaD_.25, type="l", col="steelblue",
xlab = expression(theta),
ylab = expression(p(theta~","~omega~"=.25|D")),
main = expression("Marginal"~p(theta~","~omega~"=.25|D")),
las = 1, lwd=2)
plot(theta, p_thetaD_.75, type="l", col="steelblue",
xlab = expression(theta),
ylab = expression(p(theta~","~omega~"=.75|D")),
main = expression("Marginal"~p(theta~","~omega~"=.75|D")),
las = 1, lwd=2)
}
return(invisible(list(prior = prior_matrix,
likelihood = likelihood_matrix,
posterior = posterior_matrix)))
}
'Chapter9Figure9.3' <- function(theta = seq(from=0, to=1, length=100),
omega = seq(from=0, to=1, length=100),
K = 6,
A_omega = 20,
B_omega = 20,
N = 12,
z = 9,
file_path = NULL)
{
p_ = Chapter9Figure9.2(K = K,
A_omega = A_omega,
B_omega = B_omega,
N = N,
z = z,
file_path = file_path)
return(invisible(p_))
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.