finished 3b

Oblig/3b/3b.R

@@ -0,0 +1,171 @@
+# Posterior for a Poisson process with Gamma(shape = k0, rate = t0) prior
+# and data n events observed during time t.
+poisson_posterior <- function(k0, t0, n, t) {
+  list(shape = k0 + n, rate = t0 + t)
+}
+
+# Predictive probability for future counts in a Poisson process when
+# lambda | data ~ Gamma(shape, rate).
+# R's pnbinom uses the same negative-binomial form we need here.
+predictive_count_cdf <- function(m, shape, rate, horizon) {
+  pnbinom(m, size = shape, prob = rate / (rate + horizon))
+}
+
+# Predictive probability for the k-th future waiting time.
+# We use P(T_{+k} <= l) = P(N_{+l} >= k) = 1 - P(N_{+l} <= k - 1).
+predictive_waiting_time_cdf <- function(k, l, shape, rate) {
+  1 - predictive_count_cdf(k - 1, shape, rate, l)
+}
+
+# Posterior probability P(X > Y) for independent gamma variables.
+gamma_gt_prob <- function(shape_x, rate_x, shape_y, rate_y) {
+  integrate(
+    function(x) {
+      dgamma(x, shape = shape_x, rate = rate_x) *
+        pgamma(x, shape = shape_y, rate = rate_y)
+    },
+    lower = 0,
+    upper = Inf,
+    subdivisions = 2000L,
+    rel.tol = 1e-10
+  )$value
+}
+
+# Exact probability P(U > V) for independent beta variables with integer
+# first shape parameter in U, from the finite-sum identity used in Rule A.3.3.
+beta_gt_prob_exact <- function(a, b, c, d) {
+  total <- 0
+  for (i in 0:(a - 1)) {
+    total <- total + beta(c + i, b + d) /
+      ((b + i) * beta(1 + i, b) * beta(c, d))
+  }
+  total
+}
+
+# Normal approximation for pairwise beta comparison.
+# If U and V are independent beta variables, then D = U - V is approximated
+# by a normal variable with mean E(U) - E(V) and variance Var(U) + Var(V).
+beta_compare_prob_normal <- function(a, b, c, d, direction = c("gt", "lt")) {
+  direction <- match.arg(direction)
+
+  mean_u <- a / (a + b)
+  mean_v <- c / (c + d)
+  var_u <- a * b / ((a + b)^2 * (a + b + 1))
+  var_v <- c * d / ((c + d)^2 * (c + d + 1))
+
+  mean_d <- mean_u - mean_v
+  sd_d <- sqrt(var_u + var_v)
+
+  if (direction == "gt") {
+    1 - pnorm(0, mean = mean_d, sd = sd_d)
+  } else {
+    pnorm(0, mean = mean_d, sd = sd_d)
+  }
+}
+
+# Posterior probability P(mu_x > mu_y) for Gaussian samples under neutral priors.
+# Marginally, each mean has a Student t posterior:
+# mu | data ~ t_{n-1}(xbar, s / sqrt(n)).
+gaussian_mean_gt_prob <- function(x, y) {
+  nx <- length(x)
+  ny <- length(y)
+  mx <- mean(x)
+  my <- mean(y)
+  sx <- sd(x)
+  sy <- sd(y)
+
+  integrate(
+    function(z) {
+      dt((z - mx) / (sx / sqrt(nx)), df = nx - 1) / (sx / sqrt(nx)) *
+        pt((z - my) / (sy / sqrt(ny)), df = ny - 1)
+    },
+    lower = -Inf,
+    upper = Inf,
+    subdivisions = 5000L,
+    rel.tol = 1e-10
+  )$value
+}
+
+# -----------------------------
+# Point 3a: Chapter 13, task 13d
+# -----------------------------
+
+post_13d <- poisson_posterior(k0 = 3, t0 = 73, n = 6, t = 119)
+
+cat("Chapter 13, task 13d\n")
+cat("Posterior shape =", post_13d$shape, "\n")
+cat("Posterior rate  =", post_13d$rate, "\n\n")
+
+# -----------------------------
+# Point 3b: Chapter 13, task 14c
+# -----------------------------
+
+shape_14c <- 29
+rate_14c <- 8
+k_14c <- 3
+l_14c <- 1
+m_14c <- 5
+
+p_n_le_m <- predictive_count_cdf(m = m_14c, shape = shape_14c, rate = rate_14c, horizon = l_14c)
+p_t_le_l <- predictive_waiting_time_cdf(k = k_14c, l = l_14c, shape = shape_14c, rate = rate_14c)
+
+cat("Chapter 13, task 14c\n")
+cat("P(N_{+1} <= 5) =", p_n_le_m, "\n")
+cat("P(T_{+3} <= 1) =", p_t_le_l, "\n\n")
+
+# -----------------------------
+# Point 3c: Chapter 14, task 12
+# -----------------------------
+
+odd <- c(22, 20, 21, 20, 21, 21, 19, 21)
+kjell <- c(15, 12, 32, 12, 11, 13, 14)
+
+p_mu_odd_gt_kjell <- gaussian_mean_gt_prob(odd, kjell)
+
+cat("Chapter 14, task 12\n")
+cat("Odd:   n =", length(odd), "mean =", mean(odd), "sd =", sd(odd), "\n")
+cat("Kjell: n =", length(kjell), "mean =", mean(kjell), "sd =", sd(kjell), "\n")
+cat("P(mu_Odd > mu_Kjell | data) =", p_mu_odd_gt_kjell, "\n\n")
+
+# -----------------------------
+# Point 3d: Chapter 14, task 14
+# -----------------------------
+
+p_14a <- gamma_gt_prob(7, 70, 4, 80)
+p_14b <- 1 - gamma_gt_prob(9, 20, 11, 20)
+p_14c <- 1 - gamma_gt_prob(90, 200, 110, 200)
+
+cat("Chapter 14, task 14\n")
+cat("P(Theta > Psi) in 14a =", p_14a, "\n")
+cat("P(Theta < Psi) in 14b =", p_14b, "\n")
+cat("P(Theta < Psi) in 14c =", p_14c, "\n\n")
+
+# -----------------------------
+# Point 3e: Chapter 14, task 16
+# -----------------------------
+
+# 16a: H1 is psi < pi, so exact probability is 1 - P(psi > pi).
+p_16a_exact <- 1 - beta_gt_prob_exact(2, 5, 4, 3)
+p_16a_norm <- beta_compare_prob_normal(2, 5, 4, 3, direction = "lt")
+
+# 16b: H1 is psi > pi.
+p_16b_exact <- beta_gt_prob_exact(23, 17, 17, 23)
+p_16b_norm <- beta_compare_prob_normal(23, 17, 17, 23, direction = "gt")
+
+# 16c: H1 is psi > pi.
+p_16c_exact <- beta_gt_prob_exact(20, 20, 17, 23)
+p_16c_norm <- beta_compare_prob_normal(20, 20, 17, 23, direction = "gt")
+
+# 16d: same as 16c, but all parameters multiplied by 10.
+p_16d_exact <- beta_gt_prob_exact(200, 200, 170, 230)
+p_16d_norm <- beta_compare_prob_normal(200, 200, 170, 230, direction = "gt")
+
+cat("Chapter 14, task 16\n")
+cat("16a exact  =", p_16a_exact, "\n")
+cat("16a normal =", p_16a_norm, "\n")
+cat("16b exact  =", p_16b_exact, "\n")
+cat("16b normal =", p_16b_norm, "\n")
+cat("16c exact  =", p_16c_exact, "\n")
+cat("16c normal =", p_16c_norm, "\n")
+cat("16d exact  =", p_16d_exact, "\n")
+cat("16d normal =", p_16d_norm, "\n")

Oblig/3b/oblig3b_writeup.tex

@@ -0,0 +1,486 @@
+% Paste this into your document body.
+% Preamble requirements:
+% \usepackage{amsmath, amssymb, booktabs, float, minted}
+% Compile with -shell-escape if you use minted.
+
+\section{Begreper}
+
+\subsection{Poisson-prosess}
+En Poisson-prosess er en stokastisk prosess som teller antall hendelser over tid, der
+\begin{itemize}
+  \item $N(0)=0$,
+  \item prosessen har uavhengige inkrementer,
+  \item antall hendelser i et tidsintervall av lengde $t$ er Poisson-fordelt med parameter $\lambda t$.
+\end{itemize}
+Hvis $N(t)$ er antall hendelser fram til tid $t$, sa har vi altsa
+\[
+N(t) \sim \mathrm{Poisson}(\lambda t).
+\]
+Parameteren $\lambda$ er intensiteten, det vil si forventet antall hendelser per tidsenhet.
+
+\subsection{Prediktiv fordeling}
+En prediktiv fordeling er fordelingen til en framtidig observasjon nar vi tar hensyn til
+usikkerheten i parameteren. Hvis $\theta$ er ukjent parameter og $X^+$ er en framtidig
+observasjon, er den prediktive fordelingen gitt ved
+\[
+f(x^+ \mid \text{data}) = \int f(x^+ \mid \theta) \pi(\theta \mid \text{data}) \, d\theta.
+\]
+Vi integrerer altsa observasjonsmodellen over posteriorfordelingen til parameteren.
+
+\subsection{Gamma-fordeling og gamma-gamma-fordeling}
+Gamma-fordelingen brukes ofte som prior eller posterior for en positiv parameter, for
+eksempel intensiteten $\lambda$ i en Poisson-prosess:
+\[
+\lambda \sim \Gamma(k,t).
+\]
+Hvis vi sa ser pa en framtidig ventetid eller en framtidig observasjon og integrerer ut
+$\lambda$, far vi en ny fordeling. I denne settingen far ventetiden en gamma-gamma-fordeling
+(ogsa kalt beta-prime i en alternativ parametrisering).
+
+Sammenhengen er derfor at gamma-gamma-fordelingen oppstar nar en gammafordelt parameter
+integreres ut av en gammafordelt betinget modell. Forskjellen er at gamma-fordelingen er en
+fordeling for selve parameteren, mens gamma-gamma-fordelingen er en prediktiv fordeling for
+en framtidig storrelse.
+
+\section{Oppgaver fra punkt 3}
+
+\subsection{Kapittel 13, oppgave 13d}
+Prioren er
+\[
+\lambda \sim \Gamma(k_0,t_0), \qquad k_0=3, \quad t_0=73,
+\]
+og vi observerer $n=6$ hendelser i lopet av $t=119$ tidsenheter.
+
+For en Poisson-prosess med gamma-prior er posterioren
+\[
+\lambda \mid \text{data} \sim \Gamma(k_0+n,\; t_0+t).
+\]
+Dermed far vi
+\[
+\lambda \mid \text{data} \sim \Gamma(3+6,\; 73+119)=\Gamma(9,192).
+\]
+
+\subsection{Kapittel 13, oppgave 14c}
+Vi er gitt posteriorhyperparametrene
+\[
+k_1=29, \qquad \tau_1=8,
+\]
+og skal finne
+\[
+P(T_{+3} \le 1)
+\qquad \text{og} \qquad
+P(N_{+1} \le 5).
+\]
+
+Nar $\lambda \mid \text{data} \sim \Gamma(k_1,\tau_1)$, far vi den prediktive fordelingen
+for framtidig antall hendelser som
+\[
+N_{+l} \sim \mathrm{NegBin}\!\left(k_1,\frac{\tau_1}{\tau_1+l}\right).
+\]
+Her blir det
+\[
+N_{+1} \sim \mathrm{NegBin}\!\left(29,\frac{8}{9}\right).
+\]
+Dermed er
+\[
+P(N_{+1} \le 5)
+=
+\sum_{n=0}^{5}
+\binom{n+28}{n}
+\left(\frac{8}{9}\right)^{29}
+\left(\frac{1}{9}\right)^n
+\approx 0.8298.
+\]
+
+Videre bruker vi sammenhengen
+\[
+T_{+3} \le 1
+\iff
+N_{+1} \ge 3.
+\]
+Derfor far vi
+\[
+P(T_{+3} \le 1)=P(N_{+1} \ge 3)=1-P(N_{+1}\le 2)\approx 0.6849.
+\]
+
+\subsection{Kapittel 14, oppgave 12}
+Vi bruker neutral priorer. Da far hver middelverdi en marginal posteriorfordeling av
+Student-$t$-type:
+\[
+\mu_O \mid \text{data} \sim t_{7}\!\left(\bar x_O,\frac{s_O}{\sqrt{8}}\right),
+\qquad
+\mu_K \mid \text{data} \sim t_{6}\!\left(\bar x_K,\frac{s_K}{\sqrt{7}}\right).
+\]
+
+Fra dataene far vi
+\[
+\bar x_O = 20.625, \quad s_O = 0.9161,
+\qquad
+\bar x_K = 15.5714, \quad s_K = 7.3679.
+\]
+
+Hypotesen er
+\[
+H_1:\mu_O > \mu_K,
+\]
+med signifikans $\alpha=0.2$. I denne typen oppgave velger vi $H_1$ dersom
+\[
+P(\mu_O > \mu_K \mid \text{data}) > 1-\alpha = 0.8.
+\]
+
+Ved numerisk integrasjon far vi
+\[
+P(\mu_O > \mu_K \mid \text{data}) \approx 0.9392.
+\]
+Siden $0.9392 > 0.8$, velger vi $H_1$. Dataene gir derfor tilstrekkelig stotte for at
+Odds hosteanfall i gjennomsnitt varer lenger enn Kjells.
+
+\subsection{Kapittel 14, oppgave 14}
+Vi bruker samme beslutningsregel gjennom hele oppgaven: velg $H_1$ nar den posterior
+sannsynligheten for utsagnet i $H_1$ er storre enn $1-\alpha$.
+
+\paragraph{14a}
+Vi har
+\[
+\Theta \sim \Gamma(7,70),
+\qquad
+\Psi \sim \Gamma(4,80),
+\]
+og hypotesen
+\[
+H_1:\Theta > \Psi,
+\qquad \alpha = 0.05.
+\]
+Den relevante posterior sannsynligheten er
+\[
+P(\Theta > \Psi)
+=
+\int_0^\infty f_\Theta(x)F_\Psi(x)\,dx
+\approx 0.8774.
+\]
+Siden $0.8774 < 0.95$, velger vi $H_0$.
+
+\paragraph{14b}
+Vi har
+\[
+\Theta \sim \Gamma(9,20),
+\qquad
+\Psi \sim \Gamma(11,20),
+\]
+og hypotesen
+\[
+H_1:\Theta < \Psi,
+\qquad \alpha = 0.1.
+\]
+Her gir beregningen
+\[
+P(\Theta < \Psi)\approx 0.6762.
+\]
+Siden $0.6762 < 0.9$, velger vi $H_0$.
+
+\paragraph{14c}
+Nar parameterne i 14b ganges med $10$, far vi
+\[
+\Theta \sim \Gamma(90,200),
+\qquad
+\Psi \sim \Gamma(110,200),
+\]
+med samme hypotese som i 14b. Da blir
+\[
+P(\Theta < \Psi)\approx 0.9220.
+\]
+Siden $0.9220 > 0.9$, velger vi na $H_1$.
+
+Det er altsa en forskjell: med ti ganger sa mye informasjon blir posteriorfordelingene
+mer konsentrerte, og det blir lettere a skille mellom parameterne.
+
+\subsection{Kapittel 14, oppgave 16}
+Her sammenligner vi beta-fordelte variable. Jeg viser bade eksakt beregning og
+normalapproksimasjon.
+
+For normalapproksimasjonen bruker vi at dersom
+\[
+U \sim \beta(a,b),
+\]
+sa er
+\[
+E(U)=\frac{a}{a+b},
+\qquad
+\operatorname{Var}(U)=\frac{ab}{(a+b)^2(a+b+1)}.
+\]
+For to uavhengige variable $\psi$ og $\pi$ approksimerer vi deretter
+\[
+D=\psi-\pi \approx N\!\bigl(E(\psi)-E(\pi),\operatorname{Var}(\psi)+\operatorname{Var}(\pi)\bigr).
+\]
+
+\paragraph{16a}
+\[
+\psi \sim \beta(2,5),
+\qquad
+\pi \sim \beta(4,3),
+\qquad
+H_1:\psi < \pi,
+\qquad
+\alpha=0.1.
+\]
+Eksakt beregning gir
+\[
+P(\psi < \pi) \approx 0.8788,
+\]
+mens normalapproksimasjonen gir
+\[
+P(\psi < \pi) \approx 0.8861.
+\]
+Begge er mindre enn $0.9$, sa vi velger $H_0$.
+
+\paragraph{16b}
+\[
+\psi \sim \beta(23,17),
+\qquad
+\pi \sim \beta(17,23),
+\qquad
+H_1:\psi > \pi,
+\qquad
+\alpha=0.1.
+\]
+Eksakt beregning gir
+\[
+P(\psi > \pi) \approx 0.9131,
+\]
+og normalapproksimasjonen gir
+\[
+P(\psi > \pi) \approx 0.9153.
+\]
+Begge er storre enn $0.9$, sa vi velger $H_1$.
+
+\paragraph{16c}
+\[
+\psi \sim \beta(20,20),
+\qquad
+\pi \sim \beta(17,23),
+\qquad
+H_1:\psi > \pi,
+\qquad
+\alpha=0.05.
+\]
+Eksakt beregning gir
+\[
+P(\psi > \pi) \approx 0.7520,
+\]
+og normalapproksimasjonen gir
+\[
+P(\psi > \pi) \approx 0.7527.
+\]
+Begge er mindre enn $0.95$, sa vi velger $H_0$.
+
+\paragraph{16d}
+Nar parameterne i 16c ganges med $10$, far vi
+\[
+\psi \sim \beta(200,200),
+\qquad
+\pi \sim \beta(170,230),
+\qquad
+H_1:\psi > \pi,
+\qquad
+\alpha=0.05.
+\]
+Eksakt beregning gir
+\[
+P(\psi > \pi) \approx 0.9834,
+\]
+og normalapproksimasjonen gir
+\[
+P(\psi > \pi) \approx 0.9837.
+\]
+Begge er storre enn $0.95$, sa vi velger $H_1$.
+
+Dette viser at mer data kan gi en tydelig forskjell selv nar forholdet mellom positive og
+negative observasjoner er det samme som for 16c.
+
+\subsection{R-kode}
+Listing~\ref{lst:oblig3b-r} viser R-koden som ble brukt til beregningene.
+
+\begin{listing}[H]
+\begin{minted}[fontsize=\small, linenos, breaklines, frame=lines]{r}
+# Oblig 3b calculations for points 1 and 3.
+# This script uses only ASCII characters.
+
+# -----------------------------
+# Helper functions
+# -----------------------------
+
+# Posterior for a Poisson process with Gamma(shape = k0, rate = t0) prior
+# and data n events observed during time t.
+poisson_posterior <- function(k0, t0, n, t) {
+  list(shape = k0 + n, rate = t0 + t)
+}
+
+# Predictive probability for future counts in a Poisson process when
+# lambda | data ~ Gamma(shape, rate).
+# R's pnbinom uses the same negative-binomial form we need here.
+predictive_count_cdf <- function(m, shape, rate, horizon) {
+  pnbinom(m, size = shape, prob = rate / (rate + horizon))
+}
+
+# Predictive probability for the k-th future waiting time.
+# We use P(T_{+k} <= l) = P(N_{+l} >= k) = 1 - P(N_{+l} <= k - 1).
+predictive_waiting_time_cdf <- function(k, l, shape, rate) {
+  1 - predictive_count_cdf(k - 1, shape, rate, l)
+}
+
+# Posterior probability P(X > Y) for independent gamma variables.
+gamma_gt_prob <- function(shape_x, rate_x, shape_y, rate_y) {
+  integrate(
+    function(x) {
+      dgamma(x, shape = shape_x, rate = rate_x) *
+        pgamma(x, shape = shape_y, rate = rate_y)
+    },
+    lower = 0,
+    upper = Inf,
+    subdivisions = 2000L,
+    rel.tol = 1e-10
+  )$value
+}
+
+# Exact probability P(U > V) for independent beta variables with integer
+# first shape parameter in U, from the finite-sum identity used in Rule A.3.3.
+beta_gt_prob_exact <- function(a, b, c, d) {
+  total <- 0
+  for (i in 0:(a - 1)) {
+    total <- total + beta(c + i, b + d) /
+      ((b + i) * beta(1 + i, b) * beta(c, d))
+  }
+  total
+}
+
+# Normal approximation for pairwise beta comparison.
+# If U and V are independent beta variables, then D = U - V is approximated
+# by a normal variable with mean E(U) - E(V) and variance Var(U) + Var(V).
+beta_compare_prob_normal <- function(a, b, c, d, direction = c("gt", "lt")) {
+  direction <- match.arg(direction)
+
+  mean_u <- a / (a + b)
+  mean_v <- c / (c + d)
+  var_u <- a * b / ((a + b)^2 * (a + b + 1))
+  var_v <- c * d / ((c + d)^2 * (c + d + 1))
+
+  mean_d <- mean_u - mean_v
+  sd_d <- sqrt(var_u + var_v)
+
+  if (direction == "gt") {
+    1 - pnorm(0, mean = mean_d, sd = sd_d)
+  } else {
+    pnorm(0, mean = mean_d, sd = sd_d)
+  }
+}
+
+# Posterior probability P(mu_x > mu_y) for Gaussian samples under neutral priors.
+# Marginally, each mean has a Student t posterior:
+# mu | data ~ t_{n-1}(xbar, s / sqrt(n)).
+gaussian_mean_gt_prob <- function(x, y) {
+  nx <- length(x)
+  ny <- length(y)
+  mx <- mean(x)
+  my <- mean(y)
+  sx <- sd(x)
+  sy <- sd(y)
+
+  integrate(
+    function(z) {
+      dt((z - mx) / (sx / sqrt(nx)), df = nx - 1) / (sx / sqrt(nx)) *
+        pt((z - my) / (sy / sqrt(ny)), df = ny - 1)
+    },
+    lower = -Inf,
+    upper = Inf,
+    subdivisions = 5000L,
+    rel.tol = 1e-10
+  )$value
+}
+
+# -----------------------------
+# Point 3a: Chapter 13, task 13d
+# -----------------------------
+
+post_13d <- poisson_posterior(k0 = 3, t0 = 73, n = 6, t = 119)
+
+cat("Chapter 13, task 13d\n")
+cat("Posterior shape =", post_13d$shape, "\n")
+cat("Posterior rate  =", post_13d$rate, "\n\n")
+
+# -----------------------------
+# Point 3b: Chapter 13, task 14c
+# -----------------------------
+
+shape_14c <- 29
+rate_14c <- 8
+k_14c <- 3
+l_14c <- 1
+m_14c <- 5
+
+p_n_le_m <- predictive_count_cdf(m = m_14c, shape = shape_14c, rate = rate_14c, horizon = l_14c)
+p_t_le_l <- predictive_waiting_time_cdf(k = k_14c, l = l_14c, shape = shape_14c, rate = rate_14c)
+
+cat("Chapter 13, task 14c\n")
+cat("P(N_{+1} <= 5) =", p_n_le_m, "\n")
+cat("P(T_{+3} <= 1) =", p_t_le_l, "\n\n")
+
+# -----------------------------
+# Point 3c: Chapter 14, task 12
+# -----------------------------
+
+odd <- c(22, 20, 21, 20, 21, 21, 19, 21)
+kjell <- c(15, 12, 32, 12, 11, 13, 14)
+
+p_mu_odd_gt_kjell <- gaussian_mean_gt_prob(odd, kjell)
+
+cat("Chapter 14, task 12\n")
+cat("Odd:   n =", length(odd), "mean =", mean(odd), "sd =", sd(odd), "\n")
+cat("Kjell: n =", length(kjell), "mean =", mean(kjell), "sd =", sd(kjell), "\n")
+cat("P(mu_Odd > mu_Kjell | data) =", p_mu_odd_gt_kjell, "\n\n")
+
+# -----------------------------
+# Point 3d: Chapter 14, task 14
+# -----------------------------
+
+p_14a <- gamma_gt_prob(7, 70, 4, 80)
+p_14b <- 1 - gamma_gt_prob(9, 20, 11, 20)
+p_14c <- 1 - gamma_gt_prob(90, 200, 110, 200)
+
+cat("Chapter 14, task 14\n")
+cat("P(Theta > Psi) in 14a =", p_14a, "\n")
+cat("P(Theta < Psi) in 14b =", p_14b, "\n")
+cat("P(Theta < Psi) in 14c =", p_14c, "\n\n")
+
+# -----------------------------
+# Point 3e: Chapter 14, task 16
+# -----------------------------
+
+# 16a: H1 is psi < pi, so exact probability is 1 - P(psi > pi).
+p_16a_exact <- 1 - beta_gt_prob_exact(2, 5, 4, 3)
+p_16a_norm <- beta_compare_prob_normal(2, 5, 4, 3, direction = "lt")
+
+# 16b: H1 is psi > pi.
+p_16b_exact <- beta_gt_prob_exact(23, 17, 17, 23)
+p_16b_norm <- beta_compare_prob_normal(23, 17, 17, 23, direction = "gt")
+
+# 16c: H1 is psi > pi.
+p_16c_exact <- beta_gt_prob_exact(20, 20, 17, 23)
+p_16c_norm <- beta_compare_prob_normal(20, 20, 17, 23, direction = "gt")
+
+# 16d: same as 16c, but all parameters multiplied by 10.
+p_16d_exact <- beta_gt_prob_exact(200, 200, 170, 230)
+p_16d_norm <- beta_compare_prob_normal(200, 200, 170, 230, direction = "gt")
+
+cat("Chapter 14, task 16\n")
+cat("16a exact  =", p_16a_exact, "\n")
+cat("16a normal =", p_16a_norm, "\n")
+cat("16b exact  =", p_16b_exact, "\n")
+cat("16b normal =", p_16b_norm, "\n")
+cat("16c exact  =", p_16c_exact, "\n")
+cat("16c normal =", p_16c_norm, "\n")
+cat("16d exact  =", p_16d_exact, "\n")
+cat("16d normal =", p_16d_norm, "\n")
+\end{minted}
+\caption{R-script for beregningene i Oblig 3b}
+\label{lst:oblig3b-r}
+\end{listing}