# Load packages
library(bayesrules)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
#posterior pdf of mu
Y<-c(7.1,8.9,8.4,8.6)
mean(Y)
## [1] 8.25
summarize_normal_normal(mean = 10, sd = 1.2,sigma = 1.3,y_bar =8.25, n = 4)
## model mean mode var sd
## 1 prior 10.00000 10.00000 1.4400000 1.2000000
## 2 posterior 8.64698 8.64698 0.3266577 0.5715398
#a)Utilize grid approximation with grid values μ∈{5,6,7,…,15} to approximate the posterior model of μ.
#Step 1: Define a grid of 11 mu values
grid_data <- data.frame(mu_grid = seq(from = 5, to = 15, length = 11))
#Step 2: Evaluate the prior & likelihood at each μ
grid_data <- grid_data %>%
mutate(prior = dnorm(mu_grid, 10, 1.2),
likelihood = dnorm(7.1, mu_grid) * dnorm(8.9, mu_grid) * dnorm(8.4, mu_grid) * dnorm(8.6, mu_grid))
#Step 3: Approximate the posterior
grid_data <- grid_data %>%
mutate(unnormalized = likelihood * prior,
posterior = unnormalized / sum(unnormalized))
# Examine the grid approximated posterior
round(grid_data, 2)
## mu_grid prior likelihood unnormalized posterior
## 1 5 0.00 0.00 0 0.00
## 2 6 0.00 0.00 0 0.00
## 3 7 0.01 0.00 0 0.00
## 4 8 0.08 0.01 0 0.49
## 5 9 0.23 0.00 0 0.51
## 6 10 0.33 0.00 0 0.00
## 7 11 0.23 0.00 0 0.00
## 8 12 0.08 0.00 0 0.00
## 9 13 0.01 0.00 0 0.00
## 10 14 0.00 0.00 0 0.00
## 11 15 0.00 0.00 0 0.00
# Plot the grid approximated posterior
ggplot(grid_data, aes(x = mu_grid, y = posterior)) +
geom_point() +
geom_segment(aes(x = mu_grid, xend = mu_grid, y = 0, yend = posterior))
#Set the seed
set.seed(841339)
#Step 4: sample from the discretized posterior
post_sample <- sample_n(grid_data, size = 10000,
weight = posterior, replace = TRUE)
# Histogram of the grid simulation with posterior pdf
ggplot(post_sample, aes(x = mu_grid)) +
geom_histogram(aes(y = ..density..), color = "white") +
stat_function(fun = dnorm, args = list(8.64698, 0.5715398)) +
lims(x = c(0, 15))
## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_bar()`).
#b)Repeat part a using a grid of 201 equally spaced values between 5 and 15.
#Step 1: Define a grid of 11 mu values
grid_data <- data.frame(mu_grid = seq(from = 5, to = 15, length = 201))
#Step 2: Evaluate the prior & likelihood at each μ
grid_data <- grid_data %>%
mutate(prior = dnorm(mu_grid, 10, 1.2),
likelihood = dnorm(7.1, mu_grid) * dnorm(8.9, mu_grid) * dnorm(8.4, mu_grid) * dnorm(8.6, mu_grid))
#Step 3: Approximate the posterior
grid_data <- grid_data %>%
mutate(unnormalized = likelihood * prior,
posterior = unnormalized / sum(unnormalized))
# Examine the grid approximated posterior
round(grid_data, 2)
## mu_grid prior likelihood unnormalized posterior
## 1 5.00 0.00 0.00 0 0.00
## 2 5.05 0.00 0.00 0 0.00
## 3 5.10 0.00 0.00 0 0.00
## 4 5.15 0.00 0.00 0 0.00
## 5 5.20 0.00 0.00 0 0.00
## 6 5.25 0.00 0.00 0 0.00
## 7 5.30 0.00 0.00 0 0.00
## 8 5.35 0.00 0.00 0 0.00
## 9 5.40 0.00 0.00 0 0.00
## 10 5.45 0.00 0.00 0 0.00
## 11 5.50 0.00 0.00 0 0.00
## 12 5.55 0.00 0.00 0 0.00
## 13 5.60 0.00 0.00 0 0.00
## 14 5.65 0.00 0.00 0 0.00
## 15 5.70 0.00 0.00 0 0.00
## 16 5.75 0.00 0.00 0 0.00
## 17 5.80 0.00 0.00 0 0.00
## 18 5.85 0.00 0.00 0 0.00
## 19 5.90 0.00 0.00 0 0.00
## 20 5.95 0.00 0.00 0 0.00
## 21 6.00 0.00 0.00 0 0.00
## 22 6.05 0.00 0.00 0 0.00
## 23 6.10 0.00 0.00 0 0.00
## 24 6.15 0.00 0.00 0 0.00
## 25 6.20 0.00 0.00 0 0.00
## 26 6.25 0.00 0.00 0 0.00
## 27 6.30 0.00 0.00 0 0.00
## 28 6.35 0.00 0.00 0 0.00
## 29 6.40 0.00 0.00 0 0.00
## 30 6.45 0.00 0.00 0 0.00
## 31 6.50 0.00 0.00 0 0.00
## 32 6.55 0.01 0.00 0 0.00
## 33 6.60 0.01 0.00 0 0.00
## 34 6.65 0.01 0.00 0 0.00
## 35 6.70 0.01 0.00 0 0.00
## 36 6.75 0.01 0.00 0 0.00
## 37 6.80 0.01 0.00 0 0.00
## 38 6.85 0.01 0.00 0 0.00
## 39 6.90 0.01 0.00 0 0.00
## 40 6.95 0.01 0.00 0 0.00
## 41 7.00 0.01 0.00 0 0.00
## 42 7.05 0.02 0.00 0 0.00
## 43 7.10 0.02 0.00 0 0.00
## 44 7.15 0.02 0.00 0 0.00
## 45 7.20 0.02 0.00 0 0.00
## 46 7.25 0.02 0.00 0 0.00
## 47 7.30 0.03 0.00 0 0.00
## 48 7.35 0.03 0.00 0 0.00
## 49 7.40 0.03 0.00 0 0.00
## 50 7.45 0.03 0.00 0 0.00
## 51 7.50 0.04 0.00 0 0.00
## 52 7.55 0.04 0.00 0 0.00
## 53 7.60 0.04 0.00 0 0.01
## 54 7.65 0.05 0.00 0 0.01
## 55 7.70 0.05 0.01 0 0.01
## 56 7.75 0.06 0.01 0 0.01
## 57 7.80 0.06 0.01 0 0.01
## 58 7.85 0.07 0.01 0 0.02
## 59 7.90 0.07 0.01 0 0.02
## 60 7.95 0.08 0.01 0 0.02
## 61 8.00 0.08 0.01 0 0.02
## 62 8.05 0.09 0.01 0 0.03
## 63 8.10 0.09 0.01 0 0.03
## 64 8.15 0.10 0.01 0 0.03
## 65 8.20 0.11 0.01 0 0.03
## 66 8.25 0.11 0.01 0 0.04
## 67 8.30 0.12 0.01 0 0.04
## 68 8.35 0.13 0.01 0 0.04
## 69 8.40 0.14 0.01 0 0.04
## 70 8.45 0.14 0.01 0 0.04
## 71 8.50 0.15 0.01 0 0.04
## 72 8.55 0.16 0.01 0 0.04
## 73 8.60 0.17 0.01 0 0.04
## 74 8.65 0.18 0.01 0 0.04
## 75 8.70 0.18 0.01 0 0.04
## 76 8.75 0.19 0.01 0 0.04
## 77 8.80 0.20 0.01 0 0.04
## 78 8.85 0.21 0.00 0 0.03
## 79 8.90 0.22 0.00 0 0.03
## 80 8.95 0.23 0.00 0 0.03
## 81 9.00 0.23 0.00 0 0.02
## 82 9.05 0.24 0.00 0 0.02
## 83 9.10 0.25 0.00 0 0.02
## 84 9.15 0.26 0.00 0 0.02
## 85 9.20 0.27 0.00 0 0.01
## 86 9.25 0.27 0.00 0 0.01
## 87 9.30 0.28 0.00 0 0.01
## 88 9.35 0.29 0.00 0 0.01
## 89 9.40 0.29 0.00 0 0.01
## 90 9.45 0.30 0.00 0 0.01
## 91 9.50 0.30 0.00 0 0.00
## 92 9.55 0.31 0.00 0 0.00
## 93 9.60 0.31 0.00 0 0.00
## 94 9.65 0.32 0.00 0 0.00
## 95 9.70 0.32 0.00 0 0.00
## 96 9.75 0.33 0.00 0 0.00
## 97 9.80 0.33 0.00 0 0.00
## 98 9.85 0.33 0.00 0 0.00
## 99 9.90 0.33 0.00 0 0.00
## 100 9.95 0.33 0.00 0 0.00
## 101 10.00 0.33 0.00 0 0.00
## 102 10.05 0.33 0.00 0 0.00
## 103 10.10 0.33 0.00 0 0.00
## 104 10.15 0.33 0.00 0 0.00
## 105 10.20 0.33 0.00 0 0.00
## 106 10.25 0.33 0.00 0 0.00
## 107 10.30 0.32 0.00 0 0.00
## 108 10.35 0.32 0.00 0 0.00
## 109 10.40 0.31 0.00 0 0.00
## 110 10.45 0.31 0.00 0 0.00
## 111 10.50 0.30 0.00 0 0.00
## 112 10.55 0.30 0.00 0 0.00
## 113 10.60 0.29 0.00 0 0.00
## 114 10.65 0.29 0.00 0 0.00
## 115 10.70 0.28 0.00 0 0.00
## 116 10.75 0.27 0.00 0 0.00
## 117 10.80 0.27 0.00 0 0.00
## 118 10.85 0.26 0.00 0 0.00
## 119 10.90 0.25 0.00 0 0.00
## 120 10.95 0.24 0.00 0 0.00
## 121 11.00 0.23 0.00 0 0.00
## 122 11.05 0.23 0.00 0 0.00
## 123 11.10 0.22 0.00 0 0.00
## 124 11.15 0.21 0.00 0 0.00
## 125 11.20 0.20 0.00 0 0.00
## 126 11.25 0.19 0.00 0 0.00
## 127 11.30 0.18 0.00 0 0.00
## 128 11.35 0.18 0.00 0 0.00
## 129 11.40 0.17 0.00 0 0.00
## 130 11.45 0.16 0.00 0 0.00
## 131 11.50 0.15 0.00 0 0.00
## 132 11.55 0.14 0.00 0 0.00
## 133 11.60 0.14 0.00 0 0.00
## 134 11.65 0.13 0.00 0 0.00
## 135 11.70 0.12 0.00 0 0.00
## 136 11.75 0.11 0.00 0 0.00
## 137 11.80 0.11 0.00 0 0.00
## 138 11.85 0.10 0.00 0 0.00
## 139 11.90 0.09 0.00 0 0.00
## 140 11.95 0.09 0.00 0 0.00
## 141 12.00 0.08 0.00 0 0.00
## 142 12.05 0.08 0.00 0 0.00
## 143 12.10 0.07 0.00 0 0.00
## 144 12.15 0.07 0.00 0 0.00
## 145 12.20 0.06 0.00 0 0.00
## 146 12.25 0.06 0.00 0 0.00
## 147 12.30 0.05 0.00 0 0.00
## 148 12.35 0.05 0.00 0 0.00
## 149 12.40 0.04 0.00 0 0.00
## 150 12.45 0.04 0.00 0 0.00
## 151 12.50 0.04 0.00 0 0.00
## 152 12.55 0.03 0.00 0 0.00
## 153 12.60 0.03 0.00 0 0.00
## 154 12.65 0.03 0.00 0 0.00
## 155 12.70 0.03 0.00 0 0.00
## 156 12.75 0.02 0.00 0 0.00
## 157 12.80 0.02 0.00 0 0.00
## 158 12.85 0.02 0.00 0 0.00
## 159 12.90 0.02 0.00 0 0.00
## 160 12.95 0.02 0.00 0 0.00
## 161 13.00 0.01 0.00 0 0.00
## 162 13.05 0.01 0.00 0 0.00
## 163 13.10 0.01 0.00 0 0.00
## 164 13.15 0.01 0.00 0 0.00
## 165 13.20 0.01 0.00 0 0.00
## 166 13.25 0.01 0.00 0 0.00
## 167 13.30 0.01 0.00 0 0.00
## 168 13.35 0.01 0.00 0 0.00
## 169 13.40 0.01 0.00 0 0.00
## 170 13.45 0.01 0.00 0 0.00
## 171 13.50 0.00 0.00 0 0.00
## 172 13.55 0.00 0.00 0 0.00
## 173 13.60 0.00 0.00 0 0.00
## 174 13.65 0.00 0.00 0 0.00
## 175 13.70 0.00 0.00 0 0.00
## 176 13.75 0.00 0.00 0 0.00
## 177 13.80 0.00 0.00 0 0.00
## 178 13.85 0.00 0.00 0 0.00
## 179 13.90 0.00 0.00 0 0.00
## 180 13.95 0.00 0.00 0 0.00
## 181 14.00 0.00 0.00 0 0.00
## 182 14.05 0.00 0.00 0 0.00
## 183 14.10 0.00 0.00 0 0.00
## 184 14.15 0.00 0.00 0 0.00
## 185 14.20 0.00 0.00 0 0.00
## 186 14.25 0.00 0.00 0 0.00
## 187 14.30 0.00 0.00 0 0.00
## 188 14.35 0.00 0.00 0 0.00
## 189 14.40 0.00 0.00 0 0.00
## 190 14.45 0.00 0.00 0 0.00
## 191 14.50 0.00 0.00 0 0.00
## 192 14.55 0.00 0.00 0 0.00
## 193 14.60 0.00 0.00 0 0.00
## 194 14.65 0.00 0.00 0 0.00
## 195 14.70 0.00 0.00 0 0.00
## 196 14.75 0.00 0.00 0 0.00
## 197 14.80 0.00 0.00 0 0.00
## 198 14.85 0.00 0.00 0 0.00
## 199 14.90 0.00 0.00 0 0.00
## 200 14.95 0.00 0.00 0 0.00
## 201 15.00 0.00 0.00 0 0.00
# Plot the grid approximated posterior
ggplot(grid_data, aes(x = mu_grid, y = posterior)) +
geom_point() +
geom_segment(aes(x = mu_grid, xend = mu_grid, y = 0, yend = posterior))
#Set the seed
set.seed(841339)
#Step 4: sample from the discretized posterior
post_sample <- sample_n(grid_data, size = 10000,
weight = posterior, replace = TRUE)
# Histogram of the grid simulation with posterior pdf
ggplot(post_sample, aes(x = mu_grid)) +
geom_histogram(aes(y = ..density..), color = "white") +
stat_function(fun = dnorm, args = list(8.64698, 0.5715398)) +
lims(x = c(0, 15))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_bar()`).
