-
Notifications
You must be signed in to change notification settings - Fork 0
/
Chapter2.qmd
367 lines (234 loc) · 13.2 KB
/
Chapter2.qmd
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
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
---
title: "Chapter2"
format:
html:
code-fold: false
code-tools: false
editor: source
editor_options:
chunk_output_type: console
author: Lilla Gurtner
---
```{r setup, echo = F, message = F, error=F, include = F}
# knitr global options ----
knitr::opts_chunk$set(fig.pos = 'H',
echo = T,
message = F,
warning = F,
dpi = 600,
fig.align = "center",
fig.asp = 0.62) # golden ratio
library(tidyverse)
library(tidybayes) # nice ploting
library(brms)
library(wesanderson)
theme_set(theme_tidybayes())
#wesanderson::wes_palette("AsteroidCity2")
```
# 2E1.
Which of the expressions below correspond to the statement: the probability of rain on Monday?
(1) Pr(rain)
(2) Pr(rain|Monday) *
(3) Pr(Monday|rain)
(4) Pr(rain, Monday)/ Pr(Monday) *
# 2E2.
Which of the following statements corresponds to the expression: Pr(Monday|rain)?
(1) The probability of rain on Monday.
(2) The probability of rain, given that it is Monday.
(3) The probability that it is Monday, given that it is raining. *
(4) The probability that it is Monday and that it is raining.
# 2E3.
Which of the expressions below correspond to the statement: the probability that it is Monday, given that it is raining?
(1) Pr(Monday|rain) *
(2) Pr(rain|Monday)
(3) Pr(rain|Monday) Pr(Monday)
(4) Pr(rain|Monday) Pr(Monday)/ Pr(rain) *
(5) Pr(Monday|rain) Pr(rain)/ Pr(Monday)
# 2E4.
The Bayesian statistician Bruno de Finetti (1906–1985) began his 1973 book on probability theory with the declaration: “PROBABILITY DOES NOT EXIST.” The capitals appeared in the original, so I imagine de Finetti wanted us to shout this statement. What he meant is that probability is a device for describing uncertainty from the perspective of an observer with limited knowledge; it has no objective reality. Discuss the globe tossing example from the chapter, in light of this statement. What does it mean to say “the probability of water is 0.7”?
=> Given our knowledge / ignorance about the true state of water on the planet, we think that 70% of the planet are covered in water.
# 2M1.
Recall the globe tossing model from the chapter. Compute and plot the grid approximate posterior distribution for each of the following sets of observations. In each case, assume a uniform prior for p.
(1) W, W, W
(2) W, W, W, L
(3) L, W, W, L, W, W, W
```{r 2M1, include = T}
grit_length <- 30 # steps of hypothetical water coverage on planet
n_water <- 3 # tossig result, adapt for (2) and (3)
n_total <- 3 # n tosses , adapt for (2) and (3)
# create data
d <- tibble(hypothesized_true_p_of_water = seq(from = 0, to = 1, length = grit_length),
prior = 1/grit_length) |> #this could also be 1 or any other constant
mutate(likelyhood = dbinom(n_water, size = n_total, prob = hypothesized_true_p_of_water)) |> # likelyhood of observing the "grid value" p given a prior at this location
mutate(unstand_posterior = likelyhood * prior,
stand_posterior = unstand_posterior / sum(unstand_posterior)) |>
pivot_longer(names_to = "variable", values_to = "probability", 2:5) # for plotting in the next step
d |> ggplot(aes(x = hypothesized_true_p_of_water, y = probability, color = variable)) +
geom_point() +
geom_line() +
ggtitle(paste("probability of seeing", n_water, "waters in", n_total, "world tosses.", sep = " ")) +
scale_color_manual(values = wes_palette("Darjeeling2", n = 4, type = "discrete"))
```
# 2M2.
Now assume a prior for p that is equal to zero when p \< 0.5 and is a positive constant when p ≥ 0.5. Again compute and plot the grid approximate posterior distribution for each of the sets of observations in the problem just above.
(1) W, W, W
(2) W, W, W, L
(3) L, W, W, L, W, W, W
```{r 2M2.}
grit_length <- 20# steps of hypothetical water coverage on planet
n_water <- 3 # adapt for (2) and (3)
n_total <- 3# adapt for (2) and (3)
sequence <- seq(from = 0, to = 1, length = grit_length)
grit_bigger_than_0.5 <- sum(sequence >= 0.5) #this could also be 1 or any other constant
# create data
d <- tibble(hypothesized_true_p_of_water = seq(from = 0, to = 1, length = grit_length),
prior = case_when(hypothesized_true_p_of_water < 0.5 ~ 0,
hypothesized_true_p_of_water >= 0.5 ~ 1/grit_bigger_than_0.5)) |>
mutate(likelyhood = dbinom(n_water, size = n_total, prob = hypothesized_true_p_of_water)) |> # likelyhood of observing the grid value given a prior at this location
mutate(unstand_posterior = likelyhood * prior,
stand_posterior = unstand_posterior / sum(unstand_posterior)) |>
pivot_longer(names_to = "variable", values_to = "probability", 2:5) # for plotting
d |> ggplot(aes(x = hypothesized_true_p_of_water, y = probability, color = variable)) +
geom_point() +
geom_line() +
ggtitle(paste("probability of seeing", n_water, "waters in", n_total, "world tosses.", sep = " ")) +
scale_color_manual(values = wes_palette("Darjeeling2", n = 4, type = "discrete"))
```
Note how the likelyhood does not care about the prior: it is > 0 even at proportions of water that we know are not possible according to prior knowledge.
# 2M3.
Suppose there are two globes, one for Earth and one for Mars. The Earth globe is 70% covered in water. The Mars globe is 100% land. Further suppose that one of these globes---you don't know which---was tossed in the air and produced a "land" observation. Assume that each globe was equally likely to be tossed. Show that the posterior probability that the globe was the Earth, conditional on seeing "land" (Pr(Earth\|land)), is 0.23.
*information in the text:*
If I toss earth, Pr(land) = 0.3
If I toss mars, Pr(land) = 1
Pr(tossing_earth) = Pr(tossing_mars) = 0.5
To be shown:
Pr(Earth\|land) = (Pr(land\|Earth) \* Pr(earth)) / Pr(earth and land) = 0.23
```{r 2M3.}
p_e = 0.5
p_m = 0.5
p_l_given_e = 0.3
p_l_given_m = 1
p_e_given_l <- (p_l_given_e * p_e) / (p_e * p_l_given_e + p_m * p_l_given_m)
p_e_given_l
```
# 2M4.
Suppose you have a deck with only three cards. Each card has two sides, and each side is either black or white. One card has two black sides. The second card has one black and one white side. The third card has two white sides. Now suppose all three cards are placed in a bag and shuffled. Someone reaches into the bag and pulls out a card and places it flat on a table. A black side is shown facing up, but you don't know the color of the side facing down. Show that the probability that the other side is also black is 2/3. Use the counting method (Section 2 of the chapter 2) to approach this problem. This means counting up the ways that each card could produce the observed data (a black side facing up on the table).
::: {.callout-note collapse="true"}
## counting possibilities
all-white card =\> 0
b/w card =\> 1
all-black card =\> 2
=\> total possibilities of getting the result = 3
=\> 2/3 of which are generated by the all-black card
:::
# 2M5.
Now suppose there are four cards: B/B, B/W, W/W, and another B/B. Again suppose a card is drawn from the bag and a black side appears face up. Again calculate the probability that the other side is black.
::: {.callout-note collapse="true"}
## counting possibilities
all-white card =\> 0
b/w card =\> 1
all-black card 1 =\> 2
all-black card 2 =\> 2
total: 5 ways it can be black, 4/5 have a also black backside
:::
# 2M6.
Imagine that black ink is heavy, and so cards with black sides are heavier than cards with white sides. As a result, it's less likely that a card with black sides is pulled from the bag. So again assume there are three cards: B/B, B/W, and W/W. After experimenting a number of times, you conclude that for every way to pull the B/B card from the bag, there are 2 ways to pull the B/W card and 3 ways to pull the W/W card. Again suppose that a card is pulled and a black side appears face up. Show that the probability the other side is black is now 0.5. Use the counting method, as before.
::: {.callout-note collapse="true"}
## counting probabilities
ways to get a black card:
all-white card =\> 0 \* 3 (sampling) = 0
B/W card =\> 1 \* 2 (sampling) = 2
B/B card =\> 2 \* 1 (sampling) = 2
total 4
coming from the B/B card = 2/4 = 0.5
:::
# 2M7.
Assume again the original card problem, with a single card showing a black side face up. Before looking at the other side, we draw another card from the bag and lay it face up on the table. The face that is shown on the new card is white. Show that the probability that the first card, the one showing a black side, has black on its other side is now 0.75. Use the counting method, if you can. Hint: Treat this like the sequence of globe tosses, counting all the ways to see each observation, for each possible first card.
ww bw bb
::: {.callout-note collapse="true"}
## counting probabilities
Round 1: ??
Round 2: white
to be found: p(Round 1 = black)
*different ways the data could come about:*
- if Round 1 = bw
- Round 2 = ww (side 1 & 2)
=> 1*2
- if Round 1 = bb (side 1 & 2)
- Round 2 = bw
- Round 2 = ww (side 1 & 2)
=> 2*3
=> total ways of the data = 8
=> out of these 8, 6 come from a Round 1 == bb card =>
p(round1 = bb) => 6/8 = 0.75
:::
# 2H1.
Suppose there are two species of panda bear. Both are equally common in the wild and live in the same places. They look exactly alike and eat the same food, and there is yet no genetic assay capable of telling them apart. They differ however in their family sizes. Species A gives birth to twins 10% of the time, otherwise birthing a single infant. Species B births twins 20% of the time, otherwise birthing singleton infants. Assume these numbers are known with certainty, from many years of field research. Now suppose you are managing a captive panda breeding program. You have a new female panda of unknown species, and she has just given birth to twins. What is the probability that her next birth will also be twins?
p(t\|A) = 0.1
p(t\|B) = 0.2
p(A) = p(B) = 0.5
```{r 2H1}
# probability of the panda species
p_A <- 0.5
p_B <- 0.5
p_t_given_A <- 0.1 # likelyhood of twins per species
p_t_given_B <- 0.2
# prob of having seen twins at birth 1
p_A_given_t1 <- p_t_given_A * p_A / (p_A*p_t_given_A + p_B * p_t_given_B)
p_A_new_prior <- p_A_given_t1 # new prior bc we have seen a birth, and hence the old prior p_A was updated
p_B_given_t1 <- p_t_given_B * p_B / (p_A*p_t_given_A + p_B * p_t_given_B)
p_B_new_prior <- p_B_given_t1
# probability of getting twins again at birth 2
p_t2_given_A <- p_A_new_prior * p_t_given_A
p_t2_given_B <- p_B_new_prior * p_t_given_B
result <- p_t2_given_A + p_t2_given_B # the total prob of twins at birth 2, from both possible species combined
result
```
# 2H2.
Recall all the facts from the problem above. Now compute the probability that the panda we have is from species A, assuming we have observed only the first birth and that it was twins.
```{r 2H2}
p_A_given_t1
```
# 2H3.
Continuing on from the previous problem, suppose the same panda mother has a second birth and that it is not twins, but a singleton infant. Compute the posterior probability that this panda is species A.
```{r 2H3}
p_A_new_prior # knowledge from the first birth
p_B_new_prior
p_t_given_A <- 0.1 # probab of twins at any time
p_t_given_B <- 0.2
p_nt_given_A <- 1 - p_t_given_A # probab of having a not-twin
p_nt_given_B <- 1 - p_t_given_B
p_A_given_nt2 <- (p_nt_given_A * p_A_new_prior) / (p_nt_given_A * p_A_new_prior + p_nt_given_B * p_B_new_prior)
p_A_newest_prior <- p_A_given_nt2
# for the sake of completeness
p_B_given_nt2 <- (p_nt_given_B * p_B_new_prior) / (p_nt_given_A * p_A_new_prior + p_nt_given_B * p_B_new_prior)#
p_B_newest_prior <- p_B_given_nt2
```
# 2H4.
A common boast of Bayesian statisticians is that Bayesian inference makes it easy to use all of the data, even if the data are of different types. So suppose now that a veterinarian comes along who has a new genetic test that she claims can identify the species of our mother panda. But the test, like all tests, is imperfect. This is the information you have about the test:
• The probability it correctly identifies a species A panda is 0.8.
• The probability it correctly identifies a species B panda is 0.65.
The vet administers the test to your panda and tells you that the test is positive for species A. First ignore your previous information from the births and compute the posterior probability that your panda is species A. Then redo your calculation, now using the birth data as well.
p(testA\|A) = 0.8
p(testB\|B) = 0.65
```{r 2H4a}
# without the info from the births:
p_A <- 0.5
p_B <- 0.5
p_testA_given_A <- 0.8
p_testB_given_B <- 0.65
p_A_given_testA <- (p_testA_given_A * p_A ) /
(p_testA_given_A * p_A + # correct test result
(1 - p_testB_given_B) * p_B) # false positive for A
p_A_given_testA
```
```{r 2H4b}
# now with the birth knowledge
# we can use the final state of knowledge after 2 births from 2H3:
# p_A_newest_prior
# p_B_newest_prior
p_A_given_testA_2 <- (p_testA_given_A * p_A_newest_prior ) /
(p_testA_given_A * p_A_newest_prior + # correct test result
(1-p_testB_given_B) * p_B_newest_prior) # false positive for A
p_A_given_testA_2
```