Lab 2

Analisi Bayesiana di modelli a parametro singolo

library(dplyr) # optional
library(kableExtra)

Linguaggio: R

Esempio

  1. Dati normali, media non nota, varianza nota
  2. Stima di un tasso da dati di conteggio: un esempio idealizzato

Modelli con un solo parametro (non noto)

Dati normali con media non nota (e varianza nota)

\(y=2\), \(\sigma^2=1\), \(\tau^2_0=1.5\), \(\mu_0=0,1,\ldots,5\)

norm_norm <- function(y=2, mu0=seq(0,5,1), sigmasqr=1, tau0sqr=1.5){
  
  par(mfrow=c(3,2), mar=c(2,2,1.5,0.5))

  mu <- seq(min(c(y-4*sqrt(sigmasqr), mu0-4*sqrt(tau0sqr))), 
            max(c(y+4*sqrt(sigmasqr), mu0+4*sqrt(tau0sqr))), .001)
  
  tau1sqr <- (1/sigmasqr + 1/tau0sqr)^{-1}
  mu1 <- (y/sigmasqr + mu0/tau0sqr)*tau1sqr

    for(i in 1:length(mu0)){
    plot(mu, dnorm(mu,mu0[i], tau0sqr), xlab="",ylab="", xaxs="i", yaxs="i", yaxt="n", bty="n", xlim=range(mu), ylim=c(0, .75), col=4, main=bquote(mu[0]==.(mu0[i])), type='l')
    lines(mu, dnorm(mu,y,sigmasqr), lty=1, col=2)
    lines(mu, dnorm(mu,mu1[i],tau1sqr), col="purple")
  }
}
norm_norm()

Dati di conteggio

Esempio da BDA3, pag. 45

In un dato anno, in una data città, \(y=3\) persone sono morte di asma su una popolazione di 200,000 abitanti. (quanto è una stima ‘cruda’ del tasso per 100,000 persone?) Si assuma per \(\theta\) - tasso di lungo periodo di mortalità per asma ogni 100,000 abitanti - una a priori \(Gamma(3,5)\). Il modello ossrvazionale è quindi \(Pois(2\theta)\). Calcoliamo l’a posteriori e notiamo come essa cambi rispetto all’a priori.

exposure <- function(y=3, x=2){
  
  alpha <- 3; beta <- 5
  prior_theta <- rgamma(1000, alpha, beta)
  posterior_theta <- rgamma(1000, alpha+y, beta+x)
  
  par(mfrow=c(2,1), mar=c(2,2,1.5,0.5))
  hist(prior_theta, nclass=50, xlab=expression(theta),ylab="density", main="", yaxt="n",density=-1, prob=T, cex=2, xlim=c(0,3), col=4)
  abline(v=mean(prior_theta), lwd=2)
  mtext("PRIOR", side=3)
  
  hist(posterior_theta, nclass=50, xlab=expression(theta),ylab="density", main="", yaxt="n",density=-1, prob=T, cex=2, xlim=c(0,3), col="purple")
  abline(v=mean(posterior_theta), lwd=2)
  mtext("POSTERIOR", side=3)
  
  cat(paste("rate > 1, prior:",sum(prior_theta>1)/1000, "(", round(1-pgamma(1, alpha, beta),3), ")", "posterior:", sum(posterior_theta>1)/1000, "(", round(1-pgamma(1, alpha+y, beta+x),3), ")"))
}
exposure()

rate > 1, prior: 0.134 ( 0.125 ) posterior: 0.291 ( 0.301 )