da4l_class_credible_intervals.R

library(tidyverse)
read_csv("https://raw.githubusercontent.com/agricolamz/2025_HSE_b_da4l/refs/heads/main/data/Coretta_2017_icelandic.csv") |> 
  filter(speaker == "tt01") ->
  vowels

sd_prior <- 25
sd_data <- sd(vowels$vowel.dur)
sd_post <- 1/sqrt(1/sd_prior^2 + 1/sd_data^2)
mean_prior <- 87
mean_data <- mean(vowels$vowel.dur)
mean_post <- weighted.mean(c(mean_prior, mean_data), c(1/sd_prior^2, 1/sd_data^2))

10% + 80% + 10% = 100%


cred_int_l_80 <- qnorm(0.1, mean_post, sd_post)
cred_int_h_80 <- qnorm(0.9, mean_post, sd_post)

2.5% + 95% + 2.5%  = 100%

cred_int_l_95 <- qnorm(0.025, mean_post, sd_post)
cred_int_h_95 <- qnorm(0.975, mean_post, sd_post)

vowels |>
  ggplot(aes(vowel.dur)) +
  geom_histogram(aes(y = ..density..), fill = "grey")+
  stat_function(fun = dnorm, args = list(mean_prior,  sd_prior), color = "lightblue")+
  stat_function(fun = dnorm, args = list(mean_post,  sd_post), color = "red")+
  annotate(geom = "errorbar", y = 0, xmin = cred_int_l_80, xmax = cred_int_h_80, width = 0.001, size = 2)+
  annotate(geom = "errorbar", y = -0.001, xmin = cred_int_l_95, xmax = cred_int_h_95, width = 0.001, size = 2)+
  annotate(geom = "text", y = -0.002, x = mean_post, label = "95% credible interval")+
  annotate(geom = "text", y = 0.001, x = mean_post, label = "80% credible interval")