-
Notifications
You must be signed in to change notification settings - Fork 256
/
Copy pathbayes.R
73 lines (68 loc) · 2.44 KB
/
bayes.R
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
69
70
71
72
73
#' ---
#' title: "Regression and Other Stories: Elections Economy - Bayes"
#' author: "Andrew Gelman, Jennifer Hill, Aki Vehtari"
#' date: "`r format(Sys.Date())`"
#' output:
#' html_document:
#' theme: readable
#' toc: true
#' toc_depth: 2
#' toc_float: true
#' code_download: true
#' ---
#' Demonstration of Bayesian information aggregation. See Chapter 9 in
#' Regression and Other Stories.
#'
#' -------------
#'
#+ setup, include=FALSE
knitr::opts_chunk$set(message=FALSE, error=FALSE, warning=FALSE, comment=NA)
# switch this to TRUE to save figures in separate files
savefigs <- FALSE
#' #### Load packages
library("rprojroot")
root<-has_file(".ROS-Examples-root")$make_fix_file()
#' ## Calculations
#'
#' Prior based on a previously-fitted model using economic and
#' political condition.
theta_hat_prior <- 0.524
se_prior <- 0.041
#' Survey of 400 people, of whom 190 say they will vote for the
#' Democratic candidate
n <- 400
y <- 190
#' #### Data estimate
theta_hat_data <- y/n
se_data <- sqrt((y/n)*(1-y/n)/n)
#' #### Bayes estimate
theta_hat_bayes <- (theta_hat_prior/se_prior^2 + theta_hat_data/se_data^2) / (1/se_prior^2 + 1/se_data^2)
se_bayes <- sqrt(1/(1/se_prior^2 + 1/se_data^2))
#' ## Figures
#'
#+ eval=FALSE, include=FALSE
if (savefigs) pdf(root("ElectionsEconomy/figs","prior_data_posterior_a.pdf", height=3, width=5.5))
#+
par(mar=c(3,1,1,1), mgp=c(1.5, 0.5, 0), tck=-.02)
plot(0, 0, xlim=c(0.37, 0.67), ylim=c(0, 20), xlab=expression(theta), xaxt="n", ylab="", yaxs="i", yaxt="n", bty="n", cex.lab=1.2)
axis(1, seq(0.3, 0.7, 0.1))
curve(dnorm(x, theta_hat_prior, se_prior), n=1000, add=TRUE)
text(0.588, 5, "Prior")
curve(dnorm(x, theta_hat_data, se_data), n=1000, add=TRUE)
text(0.420, 8, "Likelihood")
#+ eval=FALSE, include=FALSE
if (savefigs) dev.off()
#+ eval=FALSE, include=FALSE
if (savefigs) pdf(root("ElectionsEconomy/figs","prior_data_posterior_b.pdf", height=3, width=5.5))
#+
par(mar=c(3,1,1,1), mgp=c(1.5, 0.5, 0), tck=-.02)
plot(0, 0, xlim=c(0.37, 0.67), ylim=c(0, 20), xlab=expression(theta), xaxt="n", ylab="", yaxs="i", yaxt="n", bty="n", cex.lab=1.2)
axis(1, seq(0.3, 0.7, 0.1))
curve(dnorm(x, theta_hat_prior, se_prior), n=1000, add=TRUE, col="gray30")
text(0.588, 5, "Prior")
curve(dnorm(x, theta_hat_data, se_data), n=1000, add=TRUE, col="gray30")
text(0.420, 8, "Likelihood")
curve(dnorm(x, theta_hat_bayes, se_bayes), n=1000, add=TRUE)
text(0.525, 15, "Posterior")
#+ eval=FALSE, include=FALSE
if (savefigs) dev.off()