Lab 3 (cont)

Normale con media e varianza non note

library(dplyr)  # optional
library(ggplot2)
theme_set(theme_minimal())
library(grid)
library(gridExtra)
library(tidyr)

A priori non informativa

Distribuzione congiunta multivariata, distribuzione condizionale, distribuzione marginale, marginalizzazione e distribuzione predittiva a posteriori

Dati

y <- c(93, 112, 122, 135, 122, 150, 118, 90, 124, 114)

Statistiche sufficienti

n <- length(y)
s2 <- var(y)
my <- mean(y)

Fattorizziamo la congiunta a posteriori p(mu,sigma2|y) in p(sigma2|y)p(mu|sigma2,y) e campioniamo dalla congiunta a posteriori utilizzando questa fattorizzazione

# funzioni di supporto per campionare e valutare distribuzione chi-quadrato inversa
# scalata
rsinvchisq <- function(n, nu, s2, ...) nu * s2/rchisq(n, nu, ...)
dsinvchisq <- function(x, nu, s2) {
    exp(log(nu/2) * nu/2 - lgamma(nu/2) + log(s2)/2 * nu - log(x) * (nu/2 + 1) - (nu *
        s2/2)/x)
}

Estraiamo 1000 valori casuali da p(sigma2|y)

ns <- 1000
sigma2 <- rsinvchisq(ns, n - 1, s2)

Campioniamo da p(mu|sigma2,y)

mu <- my + sqrt(sigma2/n) * rnorm(ns)

Campioniamo dalla distribuzione predittiva p(ynew|y) per ogni estrazione di (mu, sigma)

sigma <- sqrt(sigma2)
ynew <- rnorm(ns, mu, sigma)

Per mu, sigma e ynew calcoliamo la densità su una griglia di valori
Ranges per la griglia

t1l <- c(90, 150)
t2l <- c(10, 60)
nl <- c(50, 185)
t1 <- seq(t1l[1], t1l[2], length.out = ns)
t2 <- seq(t2l[1], t2l[2], length.out = ns)
xynew <- seq(nl[1], nl[2], length.out = ns)

Calcoliamo la densità marginale esatta di mu (pm)

# multiplication by 1/sqrt(s2/n) is due to the transformation of variable
# z=(x-mean(y))/sqrt(s2/n), see BDA3 p. 21
pm <- dt((t1 - my)/sqrt(s2/n), n - 1)/sqrt(s2/n)

Stimiamo la densità marginale usando i campioni e un’approssimazione ad hoc con kernel Gaussiano (pmk)

pmk <- density(mu, adjust = 2, n = ns, from = t1l[1], to = t1l[2])$y

Facciamo la stessa cosa per sigma (densità marginale esatta (ps)

# the multiplication by 2*t2 is due to the transformation of variable z=t2^2, see
# BDA3 p. 21
ps <- dsinvchisq(t2^2, n - 1, s2) * 2 * t2

e approssimazione (psk))

psk <- density(sigma, n = ns, from = t2l[1], to = t2l[2])$y

Calcoliamo la densità predittiva esatta (p_new)

# multiplication by 1/sqrt(s2/n) is due to the transformation of variable, see BDA3
# p. 21
p_new <- dt((xynew - my)/sqrt(s2 * (1 + 1/n)), n - 1)/sqrt(s2 * (1 + 1/n))

Valutiamo la densità congiunta su una griglia (dfj$z).
Si noti che la densità che segue non è normalizzata, ma ciò è irrilevante quando si disegnano le curve di livello.

# Combine grid points into another data frame with all pairwise combinations
dfj <- data.frame(t1 = rep(t1, each = length(t2)), t2 = rep(t2, length(t1)))

dfj$z <- dsinvchisq(dfj$t2^2, n - 1, s2) * 2 * dfj$t2 * dnorm(dfj$t1, my, dfj$t2/sqrt(n))

# breaks for plotting the contours
cl <- seq(1e-05, max(dfj$z), length.out = 6)

Visualizzazione delle densità congiunta e marginali

Visualizziamo la densità congiunta e le densità marginali della distribuzione a posteriori della normale con media e varianza non note.

Creiamo un grafico della densità marginale di mu

dfm <- data.frame(t1, Exact = pm, Empirical = pmk) %>%
    gather(grp, p, -t1)

margmu <- ggplot(dfm) + geom_line(aes(t1, p, color = grp)) + coord_cartesian(xlim = t1l) +
    labs(title = "Marginal of mu", x = "", y = "") + scale_y_continuous(breaks = NULL) +
    theme(legend.background = element_blank(), legend.position = c(0.75, 0.8), legend.title = element_blank())

Creiamo un grafico della densità marginale di sigma

dfs <- data.frame(t2, Exact = ps, Empirical = psk) %>%
    gather(grp, p, -t2)

margsig <- ggplot(dfs) + geom_line(aes(t2, p, color = grp)) + coord_cartesian(xlim = t2l) +
    coord_flip() + labs(title = "Marginal of sigma", x = "", y = "") + scale_y_continuous(breaks = NULL) +
    theme(legend.background = element_blank(), legend.position = c(0.75, 0.8), legend.title = element_blank())
Coordinate system already present. Adding new coordinate system, which will replace the existing one.

Creiamo un grafico della densità congiunta

joint1labs <- c("Samples", "Exact contour")

joint1 <- ggplot() + geom_point(data = data.frame(mu, sigma), aes(mu, sigma, col = "1"),
    size = 0.1) + geom_contour(data = dfj, aes(t1, t2, z = z, col = "2"), breaks = cl) +
    coord_cartesian(xlim = t1l, ylim = t2l) + labs(title = "Joint posterior", x = "",
    y = "") + scale_y_continuous(labels = NULL) + scale_x_continuous(labels = NULL) +
    scale_color_manual(values = c("blue", "black"), labels = joint1labs) + guides(color = guide_legend(nrow = 1,
    override.aes = list(shape = c(16, NA), linetype = c(0, 1), size = c(2, 1)))) + theme(legend.background = element_blank(),
    legend.position = c(0.5, 0.9), legend.title = element_blank())

Combiniamo i grafici

# blank plot for combining the plots
bp <- grid.rect(gp = gpar(col = "white"), draw = F)

grid.arrange(joint1, margsig, margmu, bp, nrow = 2)

Visualizzazione della distribuzione fattorizzata

Visualizziamo il campionamento fattorizzato e le corrispondenti densità marginali e condizionali.

Creiamo un’altro grafico della congiunta a posteriori

# data frame for the conditional of mu and marginal of sigma
dfc <- data.frame(mu = t1, marg = rep(sigma[1], length(t1)), cond = sigma[1] + dnorm(t1,
    my, sqrt(sigma2[1]/n)) * 100) %>%
    gather(grp, p, marg, cond)
# legend labels for the following plot
joint2labs <- c("Exact contour plot", "Sample from joint\n  post.", "Cond. distribution\n of mu",
    "Sample from the\n  marg. of sigma")
joint2 <- ggplot() + geom_contour(data = dfj, aes(t1, t2, z = z, col = "1"), breaks = cl) +
    geom_point(data = data.frame(m = mu[1], s = sigma[1]), aes(m, s, color = "2")) + geom_line(data = dfc,
    aes(mu, p, color = grp), linetype = "dashed") + coord_cartesian(xlim = t1l, ylim = t2l) +
    labs(title = "Joint posterior", x = "", y = "") + scale_x_continuous(labels = NULL) +
    scale_color_manual(values = c("blue", "darkgreen", "darkgreen", "black"), labels = joint2labs) +
    guides(color = guide_legend(nrow = 2, override.aes = list(shape = c(NA, 16, NA, NA),
        linetype = c(1, 0, 1, 1)))) + theme(legend.background = element_blank(), legend.position = c(0.5,
    0.85), legend.title = element_blank())

Creiamo un’altro grafico della marginale di sigma

margsig2 <- ggplot(data = data.frame(t2, ps)) + geom_line(aes(t2, ps), color = "blue") +
    coord_cartesian(xlim = t2l) + coord_flip() + labs(title = "Marginal of sigma", x = "",
    y = "") + scale_y_continuous(labels = NULL)
Coordinate system already present. Adding new coordinate system, which will replace the existing one.

e combiniamo i grafici

grid.arrange(joint2, margsig2, ncol = 2)

Visualizzazione della distribuzione marginale di mu

Visualizziamo la distribuzione marginale di mu come una mistura di normali.

Calcoliamo le pdfs condizionali per ogni campione

condpdfs <- sapply(t1, function(x) dnorm(x, my, sqrt(sigma2/n)))
condpdfs %>%
    dim
[1] 1000 1000

Creiamo un grafico di alcune di esse

# data frame of 25 first samples
dfm25 <- data.frame(t1, t(condpdfs[1:25, ])) %>%
    gather(grp, p, -t1)

condmu <- ggplot(data = dfm25) + geom_line(aes(t1, p, group = grp), linetype = "dashed") +
    labs(title = "Conditional distribution of mu for first 25 samples", y = "", x = "") +
    scale_y_continuous(breaks = NULL)

Creiamo un grafico della loro media

dfsam <- data.frame(t1, colMeans(condpdfs), pm) %>%
    gather(grp, p, -t1)
# labels
mulabs <- c("avg of sampled conds", "exact marginal of mu")
meanmu <- ggplot(data = dfsam) + geom_line(aes(t1, p, size = grp, color = grp)) + labs(y = "",
    x = "mu") + scale_y_continuous(breaks = NULL) + scale_size_manual(values = c(2, 0.8),
    labels = mulabs) + scale_color_manual(values = c("orange", "black"), labels = mulabs) +
    theme(legend.position = c(0.8, 0.8), legend.background = element_blank(), legend.title = element_blank())

Combiniamo i grafici

grid.arrange(condmu, meanmu, ncol = 1)

Visualizzazione della distribuzione predittiva a posteriori

Visualizziamo il campionamento dalla distribuzione predittiva a posteriori.

Calcoliamo ogni pdf predittiva con mu e sigma dati

ynewdists <- sapply(xynew, function(x) dnorm(x, mu, sigma))
ynewdists %>%
    dim
[1] 1000 1000

Creiamo un grafico della congiunta a posteriori con un diagramma dalla distribuzione predittiva a posteriori, evidenziamo il primo campione, poi mostriamo tutti i campioni

# data frame of first draw from the sample of predictive pdfs along with the exact
# predictive pdf
dfp <- data.frame(xynew, draw = ynewdists[1, ], exact = p_new)
# data frame for plotting the ynew samples
dfy <- data.frame(ynew, p = 0.02 * max(ynewdists))
# legend labels
joint3labs <- c("Samples", "Exact contour")
pred1labs <- c("Sample from the predictive distribution", "Predictive distribution given the posterior sample")
pred2labs <- c("Samples from the predictive distribution", "Exact predictive distribution",
    "Predictive distribution given the posterior sample")

joint3 <- ggplot() + geom_point(data = data.frame(mu, sigma), aes(mu, sigma, col = "1"),
    size = 0.1) + geom_contour(data = dfj, aes(t1, t2, z = z, col = "2"), breaks = cl) +
    geom_point(data = data.frame(x = mu[1], y = sigma[1]), aes(x, y), color = "red") +
    coord_cartesian(xlim = t1l, ylim = t2l) + labs(title = "Joint posterior", x = "mu",
    y = "sigma") + scale_color_manual(values = c("blue", "black"), labels = joint3labs) +
    guides(color = guide_legend(nrow = 1, override.aes = list(shape = c(16, NA), linetype = c(0,
        1), size = c(2, 1)))) + theme(legend.background = element_blank(), legend.position = c(0.5,
    0.9), legend.title = element_blank())

Creiamo un grafico della prima distribuzione predittiva campionata ed il rispettivo campione ynew

pred1 <- ggplot() + geom_line(data = dfp, aes(xynew, draw, color = "2")) + geom_point(data = dfy,
    aes(ynew[1], p, color = "1")) + coord_cartesian(xlim = nl, ylim = c(0, 0.04)) + labs(title = "Posterior predictive distribution",
    x = "ytilde", y = "") + scale_y_continuous(breaks = NULL) + scale_color_manual(values = c("red",
    "blue"), labels = pred1labs) + guides(color = guide_legend(nrow = 2, override.aes = list(linetype = c(0,
    1), shape = c(16, NA), labels = pred1labs))) + theme(legend.background = element_blank(),
    legend.position = c(0.5, 0.9), legend.title = element_blank())

Creiamo un grafico con tutti gli ynews

pred2 <- ggplot() + geom_line(data = dfp, aes(xynew, draw, color = "3"), alpha = 0.5) +
    geom_line(data = dfp, aes(xynew, exact, color = "2")) + geom_point(data = dfy, aes(ynew,
    p, color = "1"), alpha = 0.05) + coord_cartesian(xlim = nl, ylim = c(0, 0.04)) + labs(x = "ytilde",
    y = "") + scale_y_continuous(breaks = NULL) + scale_color_manual(values = c("red",
    "black", "blue"), labels = pred2labs) + guides(color = guide_legend(nrow = 3, override.aes = list(linetype = c(0,
    1, 1), shape = c(16, NA, NA), alpha = c(1, 1, 0.5), labels = pred2labs))) + theme(legend.background = element_blank(),
    legend.position = c(0.5, 0.9), legend.title = element_blank())

Combiniamo i grafici

grid.arrange(joint3, pred1, bp, pred2, nrow = 2)