aula de hoje

  1. comparação de modelos: Mr. A & Mr. B
  • 1.1 Mr. A & Mr. B e o problema gaussiano
  • 1.2 as moedas de Mr. A & Mister B
  • 1.3 métodos aproximados: AIC e BIC
  1. o eclipse solar de 29 de março de 1919
  1. comparação de modelos: reta, parábola, polinômio de terceiro grau?
  • 3.1 dados e modelos
  • 3.2 seleção aproximada de modelos
  • 3.3 seleção bayesiana de modelos

Exercício

1. comparação de modelos: Mr. A & Mr. B

1.1 Mr. A & Mr. B e o problema gaussiano

MR. A & Mr. B
MR. A & Mr. B
  • Mr. A & Mr. B têm teorias sobre um modelo com um único parâmetro, \(w\):
  • Mr. A: o parâmetro \(w\) é fixo e igual a zero (prior de A)
  • Mr. B: o parâmetro \(w\) obedece a uma distribuição gaussiana de média 0 e desvio padrão \(\Sigma\) (prior de B)
  • fazemos uma medida de \(w\) e concluímos que ela é descrita por uma gaussiana de média \(\mu\) e variância \(\sigma\) (dados)
  • qual teoria devemos preferir?

comparação bayesiana dos modelos: \[ {P(A|D) \over P(B|D)} = {P(D|A) \over P(D|B)} {P(A) \over P(B)}=\]

se não temos preferência a priori por essas teorias, os modelos podem ser comparados por suas evidências:
fator de Bayes: \[ B_{AB} =\frac{P(D|A)}{P(D|B)}= \] \[ = \sqrt{1+1/z^2} \exp\Bigg[-\frac{\lambda^2}{2(1+z^2)}\Bigg] \] onde \(\lambda = \mu/\sigma\) e \(z= \sigma/\Sigma\)

Se \(B_{AB} > 1\), A é preferido, se não, B é preferido: prefere-se o modelo que tiver maior evidência

Vamos ver como o fator de Bayes varia em função desses parâmetros:

# reprodutibilidade
set.seed(123)

fator.bayes = function(lambda,z){return(sqrt(1+1/z^2)*exp(-lambda^2/2/(1+z^2)))}

# crio uma sequencia com vários valores de z
z=seq(0,10,0.001)

# valores de lambda
lambda = c(1,2,3)
cores = c('black','blue','green')
plot(1,main='fator de Bayes',type='n', xlab=expression(paste("z=",sigma,"/",Sigma)), ylab='B_AB',xlim=c(0,10),ylim=c(0,6))
for(i in 1:3){
  b=fator.bayes(lambda[i],z)
lines(z,b, col=cores[i], lwd=2)
}
abline(h = 1, col="red" )
abline(v = 1, col="salmon",lty='dotted')
text(5,5,expression(paste(lambda,'=1')), col="black" )
text(5,4,expression(paste(lambda,'=2')), col="blue" )
text(5,3,expression(paste(lambda,'=3')), col="green" )

Notem que, se z = 1, para \(\lambda\) = 2 e 3 a escolha é B, enquanto que para \(\lambda\) = 1 a escolha é A!

Vamos agora visualizar esta solução no plano (z,\(\lambda\))

inform=seq(-1,4,0.01) # log de sigma/Sigma
lambda=seq(0.01,5,0.01)
nx=length(inform)
ny=length(lambda)
b=matrix(nrow=nx,ncol=ny)
for (i in 1:nx){
    for (j in 1:ny){
    b[i,j]=fator.bayes(lambda[j],exp(inform[i]))
    }
}
image(inform,lambda,b,xlab=expression(paste("log(",sigma,"/",Sigma,')')),ylab=expression(lambda) ,main='')

# relação entre lambda e log(z) para B=1
lc = function(x){return(sqrt((1+x^2)*log(1+1/x^2)))}
ll=rep(0,nx)
for (i in 1:nx){ll[i]=lc(exp(inform[i]))}
lines(inform,ll)
text(3,0.8,'B=1', col="blue" )

1.2 as moedas de Mr. A & Mr. B

Joga-se uma moeda n=200 vezes e obtém-se c=115 caras

  • Mr. A: probabilidade de jogando a moeda sair cara é \(w=0.5\)
  • Mr. B: probabilidade de jogando a moeda sair cara é \(w\), com \(w\) uniformemente distribuído entre 0 e 1

Nesse caso as evidências são: \[P(D|A) = {\rm Binomial}(n,c,w=0.5) =\binom{n}{c}0.5^{n}\] \[ P(D|B) = \int P(D,w|B)dw = \binom{n}{c} {\rm Beta}(c+1,n-c+1) \] função beta: \[ {\rm Beta}(\alpha, \beta) = \frac{\Gamma (\alpha) \Gamma (\beta)}{\Gamma (\alpha + \beta)} \]

#dados
n = 200
caras = 115

# evidência de A
pdA = choose(n,caras)*0.5^n
pdA
## [1] 0.005955892
# parênteses: e se supomos w=115/200?
w = 115/200
# nesse caso a evidência seria
pdAl = choose(n,caras)*w^caras*(1-w)^(n-caras)
pdAl
## [1] 0.05699112
# evidência de B
pdB = choose(n,caras)*beta(caras + 1, n - caras + 1)
pdB
## [1] 0.004975124
# fator de Bayes
B_AB = pdA/pdB
B_AB
## [1] 1.197134

Esse resultado aponta uma ligeira preferência para o modelo A, mas não dá para decidir

Vamos ver como o fator de bayes varia com \(x=\) número de caras:

#dados
n = 200

# evidência de A
pdA = function(caras) choose(n,caras)*0.5^n

# evidência de B
pdB = function(caras) choose(n,caras)*beta(caras + 1, n - caras + 1)

# fator de Bayes
B_AB = function(caras) pdA(caras)/pdB(caras)

# visualização
c = seq(100,200,1)
plot(c,B_AB(c),type='l',col='salmon',lwd=3,xlab='caras',ylab='B_AB')

1.3 métodos aproximados: AIC e BIC

AIC e BIC:

Vamos supor que ajustamos a \(n\) dados um modelo que tem \(k\) parâmetros livres e máxima verossimilhança \({\hat L}\)

AIC: Akaike Information Criterion
BIC: Bayesian Information Criterion \[ AIC = 2 k - 2 \ln {\hat L} \] \[ BIC= \ln(n) k - 2 \ln {\hat L} \] o primeiro termo penaliza a complexidade do modelo (número de parâmetros) e o segundo a qualidade do ajuste: modelos com menores \(AIC\) ou \(BIC\) devem ser preferidos

Vamos ilustrar estes dois critérios com as moedas de Mr. A & Mister B

verossimilhança binomial: probabilidade de se obter \(c\) caras em \(n\) jogadas de uma moeda, onde a probabilidade de sair cara em uma jogada é \(w\): \[P(D|A) = {\rm Binomial}(n,c,w) = \binom{n}{c}w^x(1-w)^{n-c}\]

#dados
n = 200
caras = 115

# função log(verossimilhança
lne = function(num,ncaras,w){log(choose(num,ncaras)*w^ncaras*w^(num-ncaras))}

# modelo A: w=0.5
k = 0
w = 0.5
choose(n,caras)*w^caras*w^(n-caras)
## [1] 0.005955892
lnl = lne(n,caras,w)
AIC_A = 2*k-2*lnl
BIC_A = log(n)*k-2*lnl
c(lnl,AIC_A,BIC_A)
## [1] -5.123374 10.246749 10.246749
# modelo B: w=caras/n
k = 1
w=caras/n
choose(n,caras)*w^caras*w^(n-caras)
## [1] 8213247280
lnl = lne(n,caras,w)
AIC_B = 2*k-2*lnl
BIC_B = log(n)*k-2*lnl
c(lnl,AIC_B,BIC_B)
## [1]  22.82901 -43.65803 -40.35971
AIC_A - AIC_B
## [1] 53.90478
BIC_A - BIC_B
## [1] 50.60646

AIC e BIC favorecem B

2. o eclipse solar de 29 de março de 1919

o eclipse de Sobral, 1919
o eclipse de Sobral, 1919

exercício de Trotta (arXiv:1701.01467)

Durante o eclipse solar de 29 de março de 1919 duas expedições inglesas mediram a deflexão da luz das estrelas pelo Sol. A teoria da relatividade geral de Einstein previa que o ângulo de deflexão seria \[\alpha = \frac{4GM}{c^2 R}\] enquanto a teoria newtoniana previa metade deste valor.

Nessa expressão \(G\) é a constante da gravitação universal, \(c\) a velocidade da luz, \(M\) a massa do Sol e \(R\) o raio do Sol. Para Einstein \(\alpha = 1.74\) segundos de arco.

A equipe liderada por Eddington, na Ilha do Príncipe, na África, reportou \(1.61 \pm 0.40\) seg. de arco, enquanto que a liderada por Crommelin, em Sobral, no Ceará, reportou \(1.98 \pm 0.16\) seg. de arco.

Vamos usar a comparação bayesiana de hipóteses para saber se os dados favorecem Newton ou Einstein e qual observação foi a mais decisiva.

Dados:

# valores de alfa teóricos:
alfaN = 0.87 # Newton
alfaE = 1.74 # Einstein

# valores observados: média mu e desvio padrão sigma
# Eddington: 
mu_edd = 1.61
sig_edd = 0.40

# Crommelin:
mu_cro = 1.98
sig_cro = 0.16
    1. Priores de Newton e o de Einstein para \(\alpha\): vamos supor que os priores são funções delta de Dirac: \[ P(N) \sim \delta(\alpha_N), ~~~~~~~~~~P(E) \sim \delta(\alpha_E)\]
    1. se as observações podem ser descritas por gaussianas, a verossimilhança dos dados é \[P(D|w,M) \sim N(\mu,\sigma)\] para \(M = N~ou~E\)
    1. evidência dos modelos: \[P(D|M) = \int P(D|w,M)P(M) d\alpha = {1 \over \sqrt{2 \pi} \sigma} \exp{-{(\alpha_M-\mu)^2 \over 2 \sigma^2}}\] ou \[log P(D|M) \propto -log \sigma - {(\alpha_M-\mu)^2 \over 2 \sigma^2}\]
    1. fator de Bayes entre Einstein e Newton: \[log B = - {(\alpha_E-\mu)^2 \over 2 \sigma^2} + {(\alpha_N-\mu)^2 \over 2 \sigma^2}\]
    1. para as observações de Eddington e de Crommelin:
# Eddington: 
logB = -(alfaE-mu_edd)^2/(2*sig_edd^2) + (alfaN-mu_edd)^2/(2*sig_edd^2)
logB
## [1] 1.658438
# Crommelin:
logB = -(alfaE-mu_cro)^2/(2*sig_cro^2) + (alfaN-mu_cro)^2/(2*sig_cro^2)
logB
## [1] 22.93945

Eddington: \(log B = 1.658438\)

Crommelin: \(log B = 22.93945\)

Ambas as observações favorecem Einstein, mas as de Crommelin, feitas em Sobral, têm evidência maior.

3. comparação de modelos: reta, parábola ou polinômio de terceiro grau?

Vamos ajustar dois modelos diferentes a um conjunto de dados e compará-los.

3.1 os dados e os modelos

Considere os dados:

x = c(0.42,  0.72,  0.  ,  0.3 ,  0.15, 0.09,  0.19,  0.35,  0.4 ,  0.54, 0.42,  0.69, 0.2 ,  0.88,  0.03, 0.67,  0.42,  0.56,  0.14,  0.2)
y = c(0.33,  0.41, -0.22,  0.01, -0.05, -0.05, -0.12,  0.26,  0.29,  0.39,  0.31,  0.42, -0.01,  0.58, -0.2 , 0.52,  0.15,  0.32, -0.13, -0.09)

# erros
sigma = rep(0.1,length(x))

# número de pontos
N = length(x)
N
## [1] 20
# visualização
plot.com.barras.de.erro <- function(x, y, err, ylim=NULL, ...) {
  if (is.null(ylim))
    ylim <- c(min(y-err), max(y+err))
  plot(x, y, ylim=ylim, pch=19, col = 'coral')
  arrows(x, y-err, x, y+err, length=0.05, angle=90, code=3,lwd=2)
}
plot.com.barras.de.erro(x,y,sigma)

o que ajusta melhor estes dados: um modelo linear (M1), um quadrático (M2) ou um polinômio de terceiro grau na variável x?

vamos usar lm() para ajustar os modelos, calculando os coeficientes que maximizam a verossimilhança:

# ajuste linear
M1 = lm(y~x,weights = 1/sigma^2)
w1 = M1$coefficients

# ajuste quadrático
M2 = lm(y~x+I(x^2),weights = 1/sigma^2)
w2 = M2$coefficients

# ajuste de um polinômio de terceiro grau
M3 = lm(y~x+I(x^2)+I(x^3),weights = 1/sigma^2)
w3 = M3$coefficients

vamos plotar estes dois ajustes:

# vamos começar escrevendo  uma função para calcular um polinomio de coeficientes w e grau n:
f.polinomio = function(w,x){
s=w[1]
for(i in 2:length(w)){s = s+w[i]*x^(i-1)}
return(s)
}
# plotando os dados
plot.com.barras.de.erro(x,y,sigma)
#plotando os modelos
xp=seq(0,1,0.1)
y1=f.polinomio(w1,xp)
lines(xp,y1,col='blue',lwd=3)
y2=f.polinomio(w2,xp)
lines(xp,y2,col='green',lwd=3)
y3=f.polinomio(w3,xp)
lines(xp,y3,col='purple',lwd=3)
legend('topleft', legend= c("M1", "M2",'M3'),col=c("blue", "green", 'purple'), lty=1, cex=0.8)

vamos comparar as verossimilhanças; para um polinômio: \[ {\cal L(w)} = \frac{1}{(2 \pi \sigma^2)^{N/2}} \exp \Big[ \sum_{i=1}^N \frac{\Big(y_i - (w_0+w_1 x_i + w_2 x_i^2 +...)\Big)^2}{\sigma_i^2} \Big] \]

# log da verossimilhança
logL = function(w){
s = 0
for(i in 1:N) s = s + (y[i]-f.polinomio(w,x[i]))^2/sigma[i]^2
return(-s/2-N/2*log(2*pi*sigma[1]^2)) 
}

logL(w1)
## (Intercept) 
##    22.01087
logL(w2)
## (Intercept) 
##    22.94151

M2 tem um logL maior que M1, mas isso não implica que ele seja o melhor modelo!

um polinômio de grau 3, 4 … vai ter uma verossimilhança ainda maior, simplesmente porque os modelos têm mais graus de liberdade!

vamos verificar isso vendo como logL varia com o grau do polinômio:

rlogL = rep(0,9)
for(i in 1:9){
w = as.vector(lm(y~poly(x,i, raw=TRUE))$coefficients)
rlogL[i] = logL(w)
}
plot(seq(1,9),rlogL,col='blue',xlab='grau do polinômio',ylab='log(L)',type='l',lwd=3)

vamos plotar o resultado do ajuste para n= 1,2,3 e 9:

# parâmetros para um polinômio de grau 9:
w9 = as.vector(lm(y~poly(x,9, raw=TRUE))$coefficients)
# plotando os dados
plot.com.barras.de.erro(x,y,sigma)
#plot(x,y,pch=20,col='coral')
# plotando os modelos
xp=seq(0,1,0.001)
y1=f.polinomio(w1,xp)   
lines(xp,y1,col='blue',lwd=3)
y2=f.polinomio(w2,xp)
lines(xp,y2,col='green',lwd=3)
y3=f.polinomio(w3,xp)
lines(xp,y3,col='purple',lwd=3)
y9=f.polinomio(w9,xp)
lines(xp,y9,col='gold',lwd=3)
legend('topleft', legend= c("M1", "M2", 'M3', "M9"),col=c("blue", "green", 'purple', 'gold'), lty=1, cex=0.8)

vemos que para n=9 ocorre “over-fitting”: o modelo está ajustando parte do ruído - nesse caso, os valores dos parâmetros “explodem”:

w9
##  [1] -2.209723e-01 -1.929142e+00  1.392475e+02 -1.869853e+03  1.136453e+04
##  [6] -3.704581e+04  6.913875e+04 -7.388054e+04  4.192971e+04 -9.765595e+03

Mas qual modelo é o “melhor”?

3.2 seleção com comparação aproximada de modelos

Vamos considerar, inicialmente, um procedimento frequentista, a “correlação cruzada”: calcula-se um \(\chi^2\) deixando um ponto de fora e escolhe-se o modelo que minimiza este \(\chi^2\)

# função que faz validação cruzada "leave one out" 
# para um polinômio y~poly(x, grau)
crossvalidate <- function(x, y, grau)
{
  preds <- numeric(length(x))
    for(i in 1:length(x))
    {
# tira-se o ponto i dos dados    
        x.in <- x[-i]
        x.out <- x[i]
        y.in <- y[-i]
        y.out <- x[i]
        m <- lm(y.in ~ poly(x.in, grau) )
        new <- data.frame(x.in = seq(-3, 3, by=0.1))
# e faz-se a predição para esse ponto        
        preds[i]<- predict(m, newdata=data.frame(x.in=x.out))
    }
# retorna o erro quadrático:
  return(sum((y-preds)^2))
}

vamos ilustrar o método com polinômios de grau 1 a 9;

# validação cruzada dos modelos polinomiais, armazenando os erros quadráticos no objeto cv
a <- as.data.frame(matrix(0,nrow=9,ncol=1))
a <- data.frame(cross=numeric(9))
cv =rep(0,9)
for(i in 1:9)
{
  cv[i] <- crossvalidate(x, y, grau=i)
}

plot(seq(1,9,1),cv,log='y',type='l',col='blue',lwd=3,xlab='grau do polinômio',ylab='erro da CV')

cv
## [1]    0.1368173    0.1199108    0.2536709    0.2871873    0.5083920
## [6]   14.7282137   84.6017304  172.7680619 2135.0904316

CV prefere M2.

Vamos considerar agora outros procedimentos frequentistas: AIC e BIC

AIC = function(w){2*length(w)-2*logL(w)}
BIC = function(w){log(N)*length(w)-2*logL(w)}
AIC1 = AIC(w1)
AIC2 = AIC(w2)
AIC3 = AIC(w3)
c(AIC1,AIC2,AIC3)
## (Intercept) (Intercept) (Intercept) 
##   -40.02173   -39.88303   -38.35225

Como AIC1 < AIC2, e AIC1 < AIC3, AIC escolhe M1!

BIC1 = BIC(w1)
BIC2 = BIC(w2)
BIC3 = BIC(w3)
c(BIC1,BIC2,BIC3)
## (Intercept) (Intercept) (Intercept) 
##   -38.03027   -36.89583   -34.36932

e BIC também escolhe M1

3.3 seleção bayesiana de modelos

Vamos supor que M1 e M2 têm a mesma probabilidade a priori. Então podemos comparar os modelos por suas evidências.

posterior dos parâmetros \(w\): \[ P(w|D,M) = {P(D|w,M) P(w|M) \over P(D|M)}\]

evidência (supondo priores uniformes): \[ P(D|M) = \int P(D|w,M) P(w|M) d w \propto \int P(w|D,M) d w \]

Vamos calcular a evidência por integração do posterior sobre o espaço dos parâmetros \(w\): :

#prior uniforme
logprior = function(w){0}

#posterior (não normalizado) de w
logP = function(w){return(logL(w) + logprior(w))}

# cálculo das evidências por integração multidimensional
library('cubature')
f = function(w)exp(logP(w))

t0 = Sys.time()
E1 = adaptIntegrate(f, lowerLimit = c(-Inf, -Inf), upperLimit = c(Inf, Inf))$integral
Sys.time()-t0
## Time difference of 0.1847301 secs
t0 = Sys.time()
E2 = adaptIntegrate(f, lowerLimit = c(-Inf,-Inf,-Inf), upperLimit = c(Inf,Inf,Inf))$integral
Sys.time()-t0
## Time difference of 7.851795 secs
t0 = Sys.time()
E3 = adaptIntegrate(f, lowerLimit = c(-Inf,-Inf,-Inf,-Inf), upperLimit = c(Inf,Inf,Inf,Inf))$integral
Sys.time()-t0
## Time difference of 15.76325 mins
# (note como o tempo de processamento depende do número de parâmetros do modelo!)

c(log(E1),log(E2),log(E3))
## [1] 17.16499 16.31550 16.33286
#Fator de Bayes entre M1 e M2
E1/E2
## [1] 2.338472
#Fator de Bayes entre M1 e M3
E1/E3
## [1] 2.298211

o fator de Bayes privilegia M1

log(E1/E2)
## [1] 0.8494976

mas, pela tabela de Jeffreys, não dá para distinguir entre M1 e M2!

Exercício

  1. Simule 100 pontos \((x,y)\) tal que \(y = 2 x^3 + x^2 - 2x +5 + \epsilon\), com \(x\) uniformemente distribuído em \(-2 < x < 2\) e \(\epsilon \sim N(0,1)\). Use seu número USP como a semente dos números aleatórios. Ajuste polinômios de graus 1 a 7 a esses dados e verifique qual é o melhor modelo usando BIC e AIC.