-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathweek_3_lab.Rmd
93 lines (57 loc) · 3.59 KB
/
week_3_lab.Rmd
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
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
---
title: "Intro to Bayesian inference"
author: "Lab 3"
date: "25/01/2022"
output:
html_document:
df_print: paged
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
Hand in via GitHub by 9am Monday. For Q1-4, do calculations in R.
## Question 1
Consider the happiness example from the lecture, with 118 out of 129 women indicating they are happy. We are interested in estimating $\theta$, which is the (true) proportion of women who are happy. Calculate the MLE estimate $\hat{\theta}$ and 95% confidence interval.
```{r}
theta_hat <- 118/129
theta_hat
prop.test(118, 129, conf.level=0.95, correct = FALSE)
```
MLE = 0.9147, 95% confidence interval: (0.854,0.952)
## Question 2
Assume a Beta(1,1) prior on $\theta$. Calculate the posterior mean for $\hat{\theta}$ and 95% credible interval.
```{r}
theta_hat <- (118+1)/(129+2)
theta_hat
quantile(rbeta(10000,118+1,129-118+1),probs = c(0.025,0.975))
```
Posterior mean for $\hat{\theta}$ is 0.908, and 95% credible interval is (0.8525742, 0.9519564)
## Question 3
Now assume a Beta(10,10) prior on $\theta$. What is the interpretation of this prior? Are we assuming we know more, less or the same amount of information as the prior used in Question 2?
We assume we know more compared to Question 2. We tend to believe that the mean for $\hat{\theta}$ is close to 0.5.
## Question 4
Create a graph in ggplot which illustrates
- The likelihood (easiest option is probably to use `geom_histogram` to plot the histogram of appropriate random variables)
- The priors and posteriors in question 2 and 3 (use `stat_function` to plot these distributions)
Comment on what you observe.
```{r}
h = hist(c(rep(0,129-118),rep(1,118)), plot=FALSE)
h$density = h$counts/sum(h$counts)
plot(h,freq=FALSE,xlab='Happiness',ylab='Likelihood',main='Histogram of Likelihood')
```
```{r}
plot(c(0:100)/100,dbeta(c(0:100)/100,1,1),type='l',ylim=c(0,18))
lines(c(0:100)/100,dbeta(c(0:100)/100,10,10),col='red')
lines(c(0:100)/100,dbeta(c(0:100)/100,118+1,129-118+1),col='blue')
lines(c(0:100)/100,dbeta(c(0:100)/100,118+10,129-118+10),col='green')
legend(0.1,15,legend=c("Prior of Beta(1,1)", "Prior of Beta(10,10)","Posterior of Beta(1,1)", "Posterior of Beta(10,10)"),
col=c('black',"red", "blue",'green'), lty=1, cex=0.8)
```
Beta(10,10) prior drags its posterior towards 0.5.
## Question 5
(No R code required) A study is performed to estimate the effect of a simple training program on basketball free-throw shooting. A random sample of 100 college students is recruited into the study. Each student first shoots 100 free-throws to establish a baseline success probability. Each student then takes 50 practice shots each day for a month. At the end of that time, each student takes 100 shots for a final measurement. Let $\theta$ be the average improvement in success probability. $\theta$ is measured as the final proportion of shots made minus the initial proportion of shots made.
Given two prior distributions for $\theta$ (explaining each in a sentence):
- A noninformative prior, and
- A subjective/informative prior based on your best knowledge
Non-informative: Unif(-1,1). This prior has the same density for all possibilities (from -1 to 1).
Informative: something like beta(1,5). First, It's not very likeli that after training, students get worse, so I put no weight below zero. On the other hand, the improvement is not likely to be as high as 0.9 or something, so I put more weight on small positive values. Arguably, we can put some weight on the negative part.