Lab 3
Modelli multiparametrici: Normale e Multinomiale
library(dplyr) # optional
library(kableExtra)
library(ggplot2) # the last three libr for ex 2
library(gridExtra)
library(tidyr)
Esempi
Il modello Normale univariato con media e varianza non note e a priori non-informativa
Stima della velocità della luce (impianto teorico come sopra)
Il modello Multinomiale
Riferimento
- Congdon: ex. 2.4, 21, Survival times from carcinoma
- Gelman: sec. 3.2, p. 66, Estimating the speed of light
- Congdon: ex. 2.16, 39, Cancers in women
Linguaggio
- BUGS, Program 2.4 Carcinoma survival
R, norm2.r - R
- BUGS, Program 2.16 Female Cancers
R, mnomial.r
Soggetto
- Qual è la durata della sopravvivenza prevista per un nuovo paziente nelle stesse condizioni dei pazienti dell’indagine? E la probabilità di un tempo di sopravvivenza superiore a \(150\) settimane? E l’età incide su questi valori?
- Newcomb misurò (\(n=66\)) il tempo necessario alla luce per viaggiare dal suo laboratorio sul fiume Potomac a uno specchio alla base del Washington Monument e ritorno, una distanza totale di 7422 metri. Sulla base di questi dati, qual è una stima della velocità della luce?
- Nel 1995, il cancro al seno nel Regno Unito ha provocato il \(18\)% delle morti femminili per cancro. Qual è l’intervallo di credibilità per questo tasso in una situazione in cui abbiamo 11 tipi di cancro?
Esempio su normale con media e varianza non note e a priori non-informativa
Dati
dat <- read.table("norm2.dat", header=T)
par(mfrow=c(1,2))
for(i in 1:2) hist(dat[,i], main="", xlab=names(dat)[i])
Trasformiamo \(y\) in \(\log(y)\)
Analisi Bayesiana con a priori non-informativa
# Codice R per la generazione di campioni dalla densità a posteriori congiunta di (mu, sigma2) nel modello di dati Normale(mu, sigma2) con a priori non informativa p(mu,sigma2) proporzionale a 1/sigma2.
#
norm2.r <- function(y, m=10000)
{
n <- length(y)
# s2: sample variance
s2 <- 1/(n-1)*sum((y-mean(y))^2)
ybar <- mean(y)
# Simulation from the marginal posterior density for sigma2
chi <- rchisq(m, n-1)
sigma2post <- (n-1)*s2/chi
# Simulation from conditional posterior density for mu
# and from the (conditional) posterior predictive density for ynew.
# Simulation from a transformation of ynew: the probability that survival time exceeds 150 (log150=5.01)
mupostcond <- ynew <- g150 <- vector(length=m)
for(i in 1:m){
mupostcond[i] <- rnorm(1, ybar, sqrt(sigma2post[i]/n))
ynew[i] <- rnorm(1, mupostcond[i], sqrt(sigma2post[i]))
g150[i] <- as.numeric(ynew[i]>5.01)
}
# Alternatively: simulation from posterior predictive distribution for a future observation ynew;
# probability computation of g.150 mean value
ti <- rt(m,n-1)
ynew_d <- sqrt(1+1/n)*sqrt(s2)*ti+ybar
# Posterior probability of survival time exceeding 150
g150_exact <- 1-pt((log(150)-ybar)/(sqrt(1+1/n)*sqrt(s2)),n-1)
return(list(sigma=sigma2post,mu=mupostcond,ynew=ynew,g150=g150,ynew_d=ynew_d,g150_exact=g150_exact))
# return(list(sigma=sigma2post,mu=mupostcond,ynew=ynew,g150=g150))
}
out <- norm2.r(y)
summary(out$mu)
Min. 1st Qu. Median Mean 3rd Qu. Max.
2.475 3.687 3.860 3.861 4.029 4.936
sd(out$mu)
[1] 0.2580708
summary(out$ynew)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-3.579 3.062 3.846 3.846 4.619 9.679
sd(out$ynew)
[1] 1.199971
summary(out$ynew_d)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-1.366 3.082 3.862 3.862 4.643 11.491
summary(out$sigma)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.4115 1.0137 1.2602 1.3558 1.5812 8.0492
sd(out$sigma)
[1] 0.4957498
summary(out$g150)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000 0.0000 0.0000 0.1581 0.0000 1.0000
out$g150_exact
[1] 0.1598639
bru=20
par(mfrow=c(2,2), mar=c(2,2,2,0.5))
hist(out$mu, prob=T, br=bru, main=expression(mu))
lines(density(out$mu))
hist(out$sigma, prob=T, br=bru, main=expression(sigma^2))
lines(density(out$sigma))
hist(out$ynew, prob=T, br=bru, main=expression(y[new]))
lines(density(out$ynew))
#hist(out$ynew_d)
lines(density(out$ynew_d), col=2)
hist(out$mu, prob=T, br=bru, main=bquote(mu~","~ y[new]), xlim=c(0,8))
lines(density(out$mu))
lines(density(out$ynew_d), col=2)
Stima della velocità della luce utilizzando il modello normale
Dati
Terza serie di misurazioni di Newcomb del tempo di passaggio della luce, effettuate dal 24 luglio 1882 al 5 settembre 1882. I valori indicati *0,001+24,8 sono le misurazioni di Newcomb, lette lungo le colonne, registrate in milionesimi di secondo. La tabella mostra tre set di dati raccolti in giorni diversi.
In Moodle i dati sono contenuti nel file light.txt
Nota: il tempo effettivo misurato da Newcomb si ottiene a partire dai valori nei dati addizionando \(24800\) nanosecondi (\(10^{−9}\) secondi) (e.g., il primo valore osservato corrisponde a \(0.000024828\) secondi).
The data are recorded as deviations from 24800 nanoseconds.
Statistiche sufficienti
Se eliminiamo i due valori anomali (che hanno valore negativo e riteniamo così solo valori positivi) le statistiche sufficienti diventano
y_pos <- y[y > 0]
range(y_pos)
[1] 16 40
n_pos <- length(y_pos) %>% print()
[1] 64
s2_pos <- var(y_pos)
my_pos <- mean(y_pos)
cat(my_pos, sqrt(s2_pos))
27.75 5.083431
Calcoliamo la densità di mu su questi punti
Calcoliamo la densità marginale esatta per mu (si usa la relazione delle densità di variabili trasformate: la densità della var su scala naturale è uguale alla densità della variabile standardizzata moltiplicata per lo Jacobiano) (Si faccia un confronto con l’es sopra)
# multiplication by 1./sqrt(s2/n) is due to the transformation of variable
# z=(x-mean(y))/sqrt(s2/n), see BDA3 p. 21
df1$pm_mu <- dt((df1$t1 - my)/sqrt(s2/n), n-1) / sqrt(s2/n)
Calcoliamo la densità marginale esatta per mu per i dati positivi
Creiamo un istogramma delle misure
p1 <- ggplot() +
geom_histogram(aes(y), binwidth = 2, fill = 'steelblue', color = 'black') +
coord_cartesian(xlim = c(-42, 42)) +
scale_y_continuous(breaks = NULL) +
labs(title = 'Newcomb\'s measurements', x = 'y')
Creiamo un grafico del modello normale
# gather the data points into key-value pairs
df1 %>% head
t1 pm_mu pm_mu_pos
1 18.00000 6.369248e-08 1.512572e-22
2 18.01602 6.681809e-08 1.643435e-22
3 18.03203 7.009541e-08 1.785760e-22
4 18.04805 7.353169e-08 1.940564e-22
5 18.06406 7.713453e-08 2.108955e-22
6 18.08008 8.091191e-08 2.292138e-22
df2 <- gather(df1, grp, p, -t1)
df2 %>% head
t1 grp p
1 18.00000 pm_mu 6.369248e-08
2 18.01602 pm_mu 6.681809e-08
3 18.03203 pm_mu 7.009541e-08
4 18.04805 pm_mu 7.353169e-08
5 18.06406 pm_mu 7.713453e-08
6 18.08008 pm_mu 8.091191e-08
# legend labels
labs2 <- c('Posterior of mu', 'Posterior of mu given y > 0', 'Modern estimate')
p2 <- ggplot(data = df2) +
geom_line(aes(t1, p, color = grp)) +
geom_vline(aes(xintercept = 33),
linetype = 'dashed') +
coord_cartesian(xlim = c(-42, 42)) +
scale_y_continuous(breaks = NULL) +
labs(title = 'Normal model', x = 'mu', y = '') +
scale_color_manual(values = c('blue', 'darkgreen')) +
guides(color="none") +
annotate("text", x=24, y=0.26, label=labs2[1], hjust="right", size=5) +
annotate("text", x=26, y=0.58, label=labs2[2], hjust="right", size=5) +
annotate("text", x=32, y=0.7, label=labs2[3], hjust="right", size=5)
Combiniamo i grafici
Come risulta la stima della velocità della luce dai dati completi e da quelli senza i due valori anomali rispetto alla stima correntemente accettata della velocità della luce (nelle date condizioni dell’esperimento) corrispondente a \(33\) (\(0.00002483302\) secondi)?
Calcolare l’intervallo centrale per mu dai due datasets.
All 66 measurements: 23.57059 28.85365
Only positive data: 26.4802 29.0198
… This reinforces the fact that posterior inferences are only as good as the model and the experiment that produced the data. (Gelman)
Per quanto concerne i dati
il valore \(33\) corrisponde all’\(86\)esimo percentile.
Per quanto riguarda il modello, vedere nel capitolo 6 del Gelman le misure test per verificare l’adeguatezza del modello.
Dopo averne verificato l’inadeguatezza, possibili alternative sono un modello normale contaminato o un modello \(t\) di Student.
Esempio su multinomiale
Dati sui tumori femminili in UK nel 1995
# Data taken from Congdon ex. 2.16, Cancers in women
x <- c(14080,12990,6440,4350,3420,3190,2600,2420,1820,1760,23610)
Consideriamo un modello Dirichlet-Multinomial con a priori uniforme
mnomial.r <- function(x, m=10000)
{
# Simulation from the posterior density for (theta_1,...,theta_k)
pi_post <- rdiric.r(m,x+1)
# Alternatively use: rdirichlet(n, alpha)
# posterior Mean, Median and 95%CI for the category-1 probability
mm <- summary(pi_post[,1])[3:4]
ci <- quantile(pi_post[,1],c(.025,.975))
# Alternatively, analytical computation of the posterior mean
# ((x[1]+1)/sum(x+1))
return(list(mm,ci))
}
# R code for generating n random samples from a Dirichlet(a)
# It uses the fact that a Dirichlet can be obtained by sampling independent gamma densities (with equal scale parameter beta, e.g., beta=1, and shape parameter alpha equal to the corresponding element of the Dirichlet a vector), and then dividing such draws by their sum.
rdiric.r <- function(n,a) {
p <- length(a)
m <- matrix(nrow=n,ncol=p)
for (i in 1:p) {
m[,i] <- rgamma(n,a[i])
}
sumvec <- m %*% rep(1,p) # apply(m, 1 sum)
m/as.vector(sumvec)
}
# Let try rdiric.r
op <- options()
options(digits=3)
rdiric.r(1,rep(1,10))
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0.0504 0.0131 0.0505 0.0649 0.226 0.0512 0.302 0.11 0.00899 0.124
rdiric.r(4,rep(1,10))
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0.09656 0.00629 0.265019 0.00599 0.1064 0.1702 0.0500 0.24636 0.0507 0.00257
[2,] 0.10848 0.33637 0.035592 0.04567 0.0981 0.0642 0.0453 0.00881 0.0422 0.21537
[3,] 0.00433 0.10030 0.000196 0.20476 0.1071 0.1610 0.1232 0.19716 0.0351 0.06688
[4,] 0.04727 0.12242 0.079743 0.02898 0.1103 0.1798 0.1749 0.03958 0.1037 0.11333
rdiric.r(4,rep(1,10)) %>% apply(1,sum)
[1] 1 1 1 1
options(op) # reset (all) initial options
options("digits")
$digits
[1] 7