Lab 1 more

Analisi Bayesiana di Dati Binomiali

library(dplyr) # optional
library(kableExtra)
library(tidyr) # gather
library(gridExtra)

Altri esercizi riguardo alla visualzzazione e all’inferenza Bayesiana.

Uso di ggplot2

Probabilità di una nascita femminile data la condizione placenta previa

(BDA3 p. 37).

library(ggplot2)
theme_set(theme_minimal())

L’a posteriori è una Beta(438,544)

# seq creates evenly spaced values
df1 <- data.frame(theta = seq(0.375, 0.525, 0.001)) 
a <- 438
b <- 544
# dbeta computes the posterior density
df1$p <- dbeta(df1$theta, a, b)

Calcoliamo anche l’intervallo centrale del 95%

# seq creates evenly spaced values from 2.5% quantile to 97.5% quantile (i.e., 95% central interval)
# qbeta computes the value for a given quantile given parameters a and b
df2 <- data.frame(theta = seq(qbeta(0.025, a, b), qbeta(0.975, a, b), length.out = 100))
# compute the posterior density
df2$p <- dbeta(df2$theta, a, b)

Si disegni l’a posteriori (Beta(438,544)) e la linea verticale per la media della popolazione (48.8%)

ggplot2

One of the key ideas behind ggplot2 is that it allows you to easily iterate, building up a complex plot a layer at a time. Each layer can come from a different dataset and have a different aesthetic mappings, making it possible to create sophisticated plots that display data from multiple sources.

ggplot(mapping = aes(theta, p)) +
  geom_line(data = df1) + 
  # Add a layer of colorized 95% posterior interval
  geom_area(data = df2, aes(fill='1')) +
  # Add the proportion of girl babies in general population
  geom_vline(xintercept = 0.488, linetype='dotted', linewidth=1.25) +
  # Decorate the plot a little
  labs(title='Uniform prior -> Posterior is Beta(438,544)',x=expression(theta), y = '') +
  scale_y_continuous(expand = c(0, 0.1), breaks = NULL) +
  scale_fill_manual(values = 'lightblue', labels = '95% posterior interval') +
  theme(legend.position = 'bottom', legend.title = element_blank())

change from ggplot2 3.4.0

Warning: Using size aesthetic for lines was deprecated in ggplot2 3.4.0. ℹ Please use linewidth instead.

Effetto di diverse a priori in un modello binomiale

Effetto dell’a priori e confronto delle a posteriori con valori dei parametri differenti per le a priori beta

Oltre a ggplot2 per la visualizzazione, tidyr è usato per lavorare sui data frame

Dati osservati: 437 femmine e 543 maschi

a <- 437
b <- 543

La densità viene valutata nei punti equamente distribuiti tra 0.375 e 0.525

df1 <- data.frame(theta = seq(0.375, 0.525, 0.001))

A posteriori con Beta(1,1), ie. a priori uniforme

df1$pu <- dbeta(df1$theta, a+1, b+1)

Tre scelte differenti per le a priori

Beta(\(0.488*2\),\((1-0.488)*2\))
Beta(\(0.488*20\),\((1-0.488)*20\))
Beta(\(0.488*200\),\((1-0.488)*200\))

n <- c(2, 20, 200) # conteggi a priori
apr <- 0.488 # prop di successo a priori

# helperf ritorna, dati un numero di osservazioni a priori, una prop di successo a priori, 
# un numero di successi e fallimenti osservati e un dataframe con i valori per theta, 
# un dataframe nuovo con valori a priori e a posteriori calcolati nei punti theta.
# 
helperf <- function(n, apr, a, b, df)
  cbind(df, pr = dbeta(df$theta, n*apr, n*(1-apr)), po = dbeta(df$theta, n*apr + a, n*(1-apr) + b), n = n)

# lapply the above function over prior counts n, then bind by row, and lastly gather results into key-value pairs.
df2 <- lapply(n, helperf, apr, a, b, df1) %>% do.call(rbind, args = .)  %>%
  gather(key = grp, value = p, -c(theta, n), factor_key = T)
df2 %>% head; df2 %>% tail
  theta n grp            p
1 0.375 2  pu 0.0008507679
2 0.376 2  pu 0.0011418738
3 0.377 2  pu 0.0015257278
4 0.378 2  pu 0.0020295190
5 0.379 2  pu 0.0026876416
6 0.380 2  pu 0.0035433725
     theta   n grp            p
1354 0.520 200  po 0.0006885263
1355 0.521 200  po 0.0005006071
1356 0.522 200  po 0.0003622430
1357 0.523 200  po 0.0002608721
1358 0.524 200  po 0.0001869727
1359 0.525 200  po 0.0001333673
# add correct labels for plotting
df2$title <- factor(paste0('alpha/(alpha+beta)=0.488, alpha+beta=',df2$n))
levels(df2$grp)
[1] "pu" "pr" "po"
levels(df2$grp) <- c('Posterior with unif prior', 'Prior', 'Posterior')
df2 %>% head
  theta n                       grp            p
1 0.375 2 Posterior with unif prior 0.0008507679
2 0.376 2 Posterior with unif prior 0.0011418738
3 0.377 2 Posterior with unif prior 0.0015257278
4 0.378 2 Posterior with unif prior 0.0020295190
5 0.379 2 Posterior with unif prior 0.0026876416
6 0.380 2 Posterior with unif prior 0.0035433725
                                   title
1 alpha/(alpha+beta)=0.488, alpha+beta=2
2 alpha/(alpha+beta)=0.488, alpha+beta=2
3 alpha/(alpha+beta)=0.488, alpha+beta=2
4 alpha/(alpha+beta)=0.488, alpha+beta=2
5 alpha/(alpha+beta)=0.488, alpha+beta=2
6 alpha/(alpha+beta)=0.488, alpha+beta=2

Grafico delle distribuzioni

ggplot(data = df2) +
  geom_line(aes(theta, p, color = grp)) +
  geom_vline(xintercept = apr, linetype = 'dotted') +
  facet_wrap(~ title, ncol = 1) +
  labs(x = '', y = '') +
  scale_y_continuous(breaks = NULL) +
  theme(legend.position = 'bottom', legend.title = element_blank())

Inferenza basata sulla simulazione

Simuliamo campioni da una Beta(438,544), disegniamo un istogramma (con quantili), e facciamo lo stesso per una variabile trasformata.

Campioniamo dalla a posteriori Beta(438,544). Otteniamo tutte le estrazioni in una volta e salviamole in un vettore ‘theta’

a <- 438
b <- 544
theta <- rbeta(10000, a, b)

Calcoliamo l’odds ratio per tutte le estrazioni

phi <- theta/(1 - theta)

Calcoliamo i quantili 2.5% e 97.5% aprossimati usando i campioni

quantiles <- c(0.025, 0.975)
thetaq <- quantile(theta, quantiles)
phiq <- quantile(phi, quantiles)

Istogrammi con 30 bins per theta e phi

# merge the data into one data frame for plotting
df1 <- data.frame(phi,theta) %>% gather()
#df1 %>% head
# merge quantiles into one data frame for plotting
df2 <- data.frame(phi = phiq, theta = thetaq) %>% gather()
#df2
ggplot() +
  geom_histogram(data = df1, aes(value), bins = 30) +
  geom_vline(data = df2, aes(xintercept = value), linetype = 'dotted', linewidth = 1.25) +
  facet_wrap(~key, ncol = 1, scales = 'free_x')  +
  labs(x = '', y = '') +
  scale_y_continuous(breaks = NULL)

Se vogliamo aggoiungere la densità

ggplot(data = df1, aes(value)) +
  geom_histogram(aes(y = after_stat(density)), fill="white", col="black", bins = 30) +
  geom_density(col = 2, fill = 2, alpha=.2) +
  geom_vline(data = df2, aes(xintercept = value), linetype = 'dotted', linewidth = 1.25) +
  facet_wrap(~key, ncol = 1, scales = 'free')  +
  labs(x = '', y = '') +
  scale_y_continuous(breaks = NULL)

Grid e inverse-cdf sampling

Campionamento con griglia e con l’inversa della funzione di ripartizione

  • Calcoliamo la distribuzione a posteriori su una griglia discreta di punti moltiplicando la verosimiglianza e un’a priori non coniugata in ogni punto e poi normalizzando sui punti.
  • Simuliamo campioni dalla distribuzione a posteriori non standard risultante utilizzando l’inversa della cdf e la griglia discreta.

A posteriori con osservazioni (437,543) e a priori uniforme (Beta(1,1))

a <- 438
b <- 544

df1 <- data.frame(theta = seq(0.1, 1, 0.001))
df1$con <- dbeta(df1$theta, a, b)

Calcoliamo la densità di un’a priori non-coniugata nei punti discreti di una griglia.

E.g., un’a priori lineare a pezzi con una certa % di probabilità al di fuori dell’intervallo (0.388, 0.588) (si veda la figura 2.4 in BDA3)

pp <- rep(1, nrow(df1))
pi <- sapply(c(0.388, 0.488, 0.588), function(pi) which(df1$theta == pi))
pi
[1] 289 389 489
pm <- 11
pp[pi[1]:pi[2]] <- seq(1, pm, length.out = length(pi[1]:pi[2])) #%>% print
pp[pi[3]:pi[2]] <- seq(1, pm, length.out = length(pi[3]:pi[2])) #%>% print

Normalizziamo la a priori

df1$nc_p <- pp / sum(pp)

Calcoliamo l’a posteriori non-normalizzata e non-coniugata sulla griglia, quindi la normalizziamo

po <- dbeta(df1$theta, a, b) * pp

df1$nc_po <- po / sum(po)

Disegniamo l’a posteriori con a priori uniforme, l’a priori non-coniugata e la corrispondente a posteriori non-coniugata

# gather the data frame into key-value pairs
# and change variable names for plotting
df2 <- gather(df1, grp, p, -theta, factor_key = T) #%>%
levels(df2$grp) <- c('Posterior with uniform prior', 'Non-conjugate prior', 'Non-conjugate posterior')
ggplot(data = df2) +
  geom_line(aes(theta, p)) +
  facet_wrap(~grp, ncol = 1, scales = 'free_y') +
  coord_cartesian(xlim = c(0.35,0.6)) +
  scale_y_continuous(breaks=NULL) +
  labs(x = '', y = '')

Inverse cdf sampling

calcoliamo la densità cumulata (funzione di ripartizione) sulla griglia

df1$cs_po <- cumsum(df1$nc_po)

Campioniamo dalla distr uniforme U(0,1)

# set.seed(seed) is used to set seed for the random number generator
set.seed(2601)
# runif(k) returns k uniform random numbers from interval [0,1]
r <- runif(10000)

Campionamento tramite l’inversa della cdf

# function to find the smallest value theta at which the cumulative sum of the posterior densities is greater than r.
invcdf <- function(r, df) df$theta[sum(df$cs_po < r) + 1]
# sapply function for each sample r. The returned values s are now random draws from the distribution.
s <- sapply(r, invcdf, df1)

Creiamo tre grafici: p1 è l’a posteriori, p2 è la cdf dell’a posteriori e p3 è l’istogramma delle estrazioni dall’a posteriori (estratte usando l’inversa-cdf)

p1 <- ggplot(data = df1) +
  geom_line(aes(theta, nc_po)) +
  coord_cartesian(xlim = c(0.35, 0.6)) +
  labs(title = 'Non-conjugate posterior', x = '', y = '') +
  scale_y_continuous(breaks = NULL)
p2 <- ggplot(data = df1) +
  geom_line(aes(theta, cs_po)) +
  coord_cartesian(xlim = c(0.35, 0.6)) +
  labs(title = 'Posterior-cdf', x = '', y = '') +
  scale_y_continuous(breaks = NULL)
p3 <- ggplot() +
  geom_histogram(aes(s, y = after_stat(density)), binwidth = 0.003) +
  geom_line(aes(s, y = after_stat(density)), stat = "bin", binwidth = 0.003, col='red') +
  coord_cartesian(xlim = c(0.35, 0.6)) +
  labs(title = 'Histogram of posterior samples', x = '', y = '') +
  scale_y_continuous(breaks = NULL)
# combine the plots
grid.arrange(p1, p2, p3)