aula 7

aula de hoje

  1. comparação de modelos: Mr. A & Mr. B
  1. a moeda é honesta?
  1. o eclipse solar de 29 de março de 1919
  1. comparação de modelos: reta, parábola, polinômio de terceiro grau?
  1. quando métodos bayesianos e frequentistas divergem
  1. Exercício
# tempo de processamento do script
t00 = Sys.time()

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 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. a moeda é honesta?

Problema: temos um conjunto de jogadas de cara/coroa e temos 3 modelos para P(cara) = \(\theta\):

  • \(H_1\): moeda honesta, \(\theta\) = 0.50
  • \(H_2\): moeda moderadamente viesada, \(\theta\) = 0.70
  • \(H_3\): moeda muito viesada, \(\theta\) = 0.85

qual modelo é o mais plausível frente aos dados?

Vamos jogar moedas, assumindo um P(cara) verdadeiro de 0.65, e atualizar as probabilidades desses três modelos após cada jogada.

Vamos usar logaritmos para calcular probabilidades. Por que?

  • probabilidades como 0,740 podem se tornar muito pequenas se multiplicadas sucessivamente, causando underflow
  • somar logaritmos de probabilidade é numericamente mais estável do que multiplicar probabilidades.

probabilidades de cada modelo: \(\theta_1=0.50,~\theta_2=0.70,~\theta_3=0.85\)

priores de cada modelo: \(P(H_1)=P(H_2)=P(H_3)=1/3.\)

verossimilhança de 1 jogada de moeda (\(x\), \(x=0\) ou 1)
verossimilhança de Bernoulli: \[P(X=x|\theta) = \theta^x (1-\theta)^{1-x}\] ou, em log, \[log P(X=x|\theta) = x log\theta+(1-x) log (1-\theta)\]

Comparação bayesiana de modelos: \[P(H_j|D) = {P(D|H_j)P(H_j) \over \sum_{m=1}^3 P(D|H_m)P(H_m)}\]

# -------------------------------
# 0. reprodutibilidade
# -------------------------------
set.seed(42)

# -------------------------------
# 1. vamos gerar dados "verdadeiros"
# -------------------------------
theta_verdadeiro <- 0.65

# número de jogadas
n_jogadas <- 40

# dados simulados:
# 1 = caras, 0 = coroas
jogadas <- rbinom(n_jogadas, size = 1, prob = theta_verdadeiro)

head(jogadas)
[1] 0 0 1 0 1 1
# vamos comparar 3 modelos:
# H1: moeda honesta            theta = 0.50
# H2: moderadamente viesada    theta = 0.70
# H3: muito viesada            theta = 0.85
prob_modelo <- c(H1 = 0.50, H2 = 0.70, H3 = 0.85)

# vamos supor que cada um desses modelos têm a mesma probabilidade a priori
prior_modelo <- c(H1 = 1/3, H2 = 1/3, H3 = 1/3)

# log da verossimilhança de uma jogada dado theta
# processo de Bernoulli: experimento aleatório que possui apenas dois 
# resultados possíveis: "sucesso" ou "fracasso"
# a probabilidade de sair cara é theta e de sair coroa é (1-theta)
log_ver_bernoulli <- function(x, theta) {
  if (x == 1) {
    log(theta)
  } else {
    log(1 - theta)
  }
}

# Bayes sequencial:
# entrada:
#   log_prior     = log dos priores dos vários modelos
#   x             = um dado (1 ou 0)
#   theta_modelos  = vetor com as probabilidades do modelo
# saída:
#   log do posterior depois de um novo dado x
update_posterior <- function(log_prior, x, theta_modelos) {
  log_post_nao_normalizado <- log_prior + sapply(theta_modelos, function(th) {
    log_ver_bernoulli(x, th)
  })

  # truque "Log-sum-exp" para estabilidade numérica
  m <- max(log_post_nao_normalizado)
  log_evidencia <- m + log(sum(exp(log_post_nao_normalizado - m)))
  log_post <- log_post_nao_normalizado - log_evidencia

  log_post
}

# atualização
log_prior <- log(prior_modelo)

posterior_historia <- matrix(
  NA_real_,
  nrow = n_jogadas + 1,
  ncol = length(prob_modelo),
  dimnames = list(NULL, names(prob_modelo))
)
posterior_historia[1, ] <- prior_modelo

for (i in seq_len(n_jogadas)) {
  log_prior <- update_posterior(log_prior, jogadas[i], prob_modelo)
  posterior_historia[i + 1, ] <- exp(log_prior)
}

# probabilidades finais do posterior:
# theta_verdadeiro
theta_verdadeiro
[1] 0.65
# Probabilidades posteriores finais após a comparação de modelos:
posterior_historia[n_jogadas + 1, ]
          H1           H2           H3 
9.703162e-01 2.968240e-02 1.375045e-06 
#visualização
matplot(
  0:n_jogadas,
  posterior_historia,
  type = "l",
  lty = 1,
  lwd = 2,
  xlab = "Número das jogadas",
  ylab = "Posterior",
  main = "comparação de modelos com Bayes sequencial",
  ylim = c(0, 1)
)
legend(
  20,0.7,
  legend = paste0(names(prob_modelo), " (theta = ", prob_modelo, ")"),
  col = seq_len(ncol(posterior_historia)),
  lty = 1,
  lwd = 2,
  bty = "n"
)
grid()

3. o eclipse solar de 29 de março de 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.

4. 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.

4.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}{2 \sigma^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

Qual é a consequência do over-fitting? O modelo não generaliza bem, isto é, sua capacidade de predição para novos valores de \(x\) que não estão nos dados fica muito ruim.

Mas qual modelo é o “melhor”? Que grau do polinômio se deve escolher?

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

Vamos considerar, inicialmente, um procedimento frequentista, a “validaçã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

4.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.1009357 secs
t0 = Sys.time()
E2 = adaptIntegrate(f, lowerLimit = c(-Inf,-Inf,-Inf), upperLimit = c(Inf,Inf,Inf))$integral
Sys.time()-t0
Time difference of 4.368353 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 9.087035 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

M1 tem maior evidência, mas, pela tabela de Jeffreys, não dá para distinguir entre M1 e M2!

5. quando os métodos bayesianos e frequentistas divergem…

Vamos discutir um exemplo similar ao considerado por Bayes em 1763.

http://jakevdp.github.io/blog/2014/06/06/frequentism-and-bayesianism-2-when-results-differ/

Carol joga bolas, de costas e sem viés, numa mesa de bilhar que tem uma marca: se elas caem de um lado da marca, Alice ganha um ponto, se caem do outro, Bob ganha um ponto; ganha o jogo quem primeiro fizer 6 pontos

num certo jogo, após 8 bolas, Alice tem 5 pontos e Bob tem 3

qual é a probabilidade de Bob ganhar o jogo?

abordagem frequentista:

  • a probabilidade \(p\) da bola cair do lado da Alice é \(p=5/8\)
  • para Bob ganhar o jogo, ele tem que marcar 3 pontos seguidos e a probabilidade disso é \[P_{freq} = (1-p)^3 = 0.0527\]

abordagem bayesiana:

  • seja \(B\) o evento ``Bob ganha’’
  • dados \(D= \{n_A,n_B\} = \{5,3\}\)
  • \(p\): probabilidade (desconhecida) de que a bola caia na área da Alice
  • queremos \(P(B|D)\)
  • note que o valor de \(p\) não interessa! ele é um nuisance parameter

usando a técnica de marginalizar sobre parâmetros que não são de interesse (no caso p): \[ P(B|D) = \int P(B,p|D)dp = \int P(B|p,D)P(p,D)dp = \frac{\int P(B|p,D) P(D|p)P(p) dp}{\int P(D|p)P(p) dp}\]

  • para ganhar a partida, Bob tem que ganhar 3 jogadas seguidas: \(P(B|p,D)=P(B|p)=(1-p)^3\)
  • vamos supor o prior \(P(p)\) uniforme entre 0 e 1
  • verossimilhança: distribuição binomial para n jogadas e k sucessos: \[P(D|p) \propto p^k (1-p)^{n-k} = p^5(1-p)^3\]

logo, \[P(B|D) = \frac{\int_0^1 (1-p)^6 p^5 dp}{\int_0^1 (1-p)^3 p^5 dp}\]

a função \(\beta\) é definida como \[\beta (n,m) = \int_0^1 (1-p)^{n-1} p^{m-1} dp\] e, portanto, \[ P(B|D) = \beta(6+1,5+1)/\beta(3+1,5+1)\simeq 0.091\]

notem que \(P_{freq} \simeq 0.053\) e \(P_{bayes} \simeq 0.091\)!

Quem tem razão? Vamos fazer uma simulação:

# tempo de processamento
t0 = Sys.time()

#reprodutibilidade
set.seed(123)

# primeiro vamos considerar os resultados analíticos:
# abordagem frequentista:
# a probabilidade p da bola cair do lado da Alice é
p=5/8
p
[1] 0.625
# para Bob ganhar o jogo, ele tem que marcar 3 pontos seguidos
# probabilidade disso
Pfreq = (1-p)^3
Pfreq
[1] 0.05273438
# abordagem bayesiana
Pbayes = beta(6+1,5+1)/beta(3+1,5+1)
Pbayes
[1] 0.09090909
# vamos agora fazer uma simulação do jogo

# dado um p (número aleatório entre 0 e 1) e uma jogada j (número aleatório entre 0 e 1): se j < p, A ganha,       se não, B ganha

# cada partida precisa de no máximo 11 jogadas para um jogador ganhar 6 vezes

# vamos determinar o número de partidas que satisfaz os dados, (A ganha, B ganha)=(5,3), e em  quantas delas B ganha

# reprodutibilidade
set.seed(123)
# n jogos com p amostrado uniformemente entre 0 e 1
n = 100000
p = runif(n)
nBganha = 0
n53 = 0

for(i in 1:n){
# modelo gerativo
jogada = runif(11)
A = cumsum(jogada < p[i])   # número de jogadas que A ganha
B = cumsum(jogada >= p[i])  # número de jogadas que B ganha

# ABC
# número de simulações iguais aos dados
if(A[8] == 5 & B[8] == 3)n53=n53+1

# número de simulações onde B ganha
if(A[8] == 5 & B[8] == 3 & max(B) >= 6)nBganha = nBganha +1
}

fB = nBganha/n53
fB
[1] 0.0917892
# tempo de processamento
Sys.time()-t0
Time difference of 0.1751797 secs

As soluções frequentista e bayesiana são diferentes, e a bayesiana é a correta!

6. Exercício

  1. Inicialize o gerador de números aleatórios com seu número USP. 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)\). Ajuste polinômios de graus 1 a 7 a esses dados e verifique qual é o melhor modelo usando BIC e AIC.
# tempo de processamento do script
Sys.time()-t00
Time difference of 9.186469 mins