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. comparação de modelos: reta, parábola, polinômio de terceiro grau?
  • 2.1 dados e modelos
  • 2.2 seleção aproximada de modelos
  • 2.3 seleção bayesiana de modelos
  1. distância de uma estrela via paralaxe
  • 3.1 distância de uma estrela
  • 3.2 a verossimilhança
  • 3.3 escolha do prior
  • 3.4 posterior
  • 3.5 exemplo prático
  • 3.6 a evidência

Exercícios

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 usar o teste p com H0: \(w = 0.5\):

# valor p
# teste p H0 = 0.5 p(>= caras)
1-pbinom(caras,n,0.5)
## [1] 0.0140627

H0 seria rejeitada com CL=0.05 mas não com CL=0.01

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

2.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”?

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

2.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.2064929 secs
t0 = Sys.time()
E2 = adaptIntegrate(f, lowerLimit = c(-Inf,-Inf,-Inf), upperLimit = c(Inf,Inf,Inf))$integral
Sys.time()-t0
## Time difference of 10.02533 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 20.75243 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!

3. Inferência de distâncias de estrelas a partir da paralaxe

O que segue é baseado no artigo Practical Guidance for Bayesian Inference in Astronomy, por Eadie et al. (arXiv:2302.04703), mas com algumas modificações e adições.

O artigo se propõe a ser um guia prático do uso da inferência bayesiana e considera a estimativa de distância de uma estrela e de um aglomerado estelar a partir de medidas de paralaxe. Vou considerar aqui o caso da determinação da distância de uma estrela.

3.1 Distância a uma estrela

A relação entre a distância \(d\) em kpc e a paralaxe \(\tilde{\omega}\) em milisegundos de arco (mas) é \[ d[kpc] = {1 \over \tilde{\omega}[mas]} \]

Nossa observação (dados) é a paralaxe da estrela e seu erro, \[ {\cal D} = \tilde{\omega} \pm \sigma, \] e queremos determinar a distância \(d\).

Para isso vamos escrever o teorema de Bayes como \[ P(d|D) = {P(D|d)P(d) \over P(D)} \]

3.2 A verossimilhança

Supondo erros gaussianos, a verossimilhança bayesiana pode ser escrita como: \[ P(D|d) = {1 \over \sqrt{2 \pi \sigma}} {\rm {\rm exp}}\Big[ -{(\tilde{\omega} -1/d)^2 \over 2 \sigma^2} \Big] \] e é uma função de distribuição de probabilidades.

A verossimilhança frequentista é escrita como função dos parâmetros (e não é uma pdf). Se consideramos que o parâmetro é \(\tilde{\omega}\), \[ {\cal L}(\tilde{\omega}) \propto {\rm {\rm exp}}\Big[ -{(\tilde{\omega} -1/d)^2 \over 2 \sigma^2} \Big]\] ou, se consideramos que o parâmetro é \(d\): \[ {\cal L}(d) \propto {\rm exp}\Big[ -{(\tilde{\omega} -1/d)^2 \over 2 \sigma^2} \Big]\]

Note que, embora essas verossimilhanças tenham a mesma expressão, em termos do parâmetro considerado elas têm formas muito diferentes. Vamos ilustrar com valores e unidades arbitrárias:

# reprodutibilidade
set.seed(123)

# dois gráficos numa linha:
par(mfrow=c(1,2))

# paralaxe 
x = seq(16,24,0.1)
L = dnorm(x,mean=20,sd=1)
# xaxt='n', yaxt='n': removem os valores nos gráficos
plot(x,L,xlab='paralaxe [mas]', ylab=expression(paste('L(',omega,')')), type='l', col='salmon', xaxt='n', yaxt='n', lwd=3)
legend(15,0.4,legend='gaussiana',bty='n')

#distância
x=seq(0.1,100,0.01)
L = exp(-(21-1000/x)^2/40)
plot(x,L,xlab='distância (pc)', ylab='L(d)', col='salmon', xaxt='n', yaxt='n', type='l', lwd=3)

3.3 Escolha do prior

Para determinarmos o posterior do parâmetro \(d\) precisamos escolher o prior. Vamos considerar algumas opções:

P1: prior não informativo em \(\tilde{\omega}\): \(P(\tilde{\omega}) = cte\) para \(\tilde{\omega} > 0\), ou \[ P(\tilde{\omega}) \sim {\rm Unif}(0, +\infty) \] onde \({\rm Unif}()\) é uma distribuição uniforme.

ou, numa versão mais geral, fracamente informativa em \(\tilde{\omega}\): \[ P(\tilde{\omega}) \sim {\rm Unif}(\tilde{\omega}_{min},\tilde{\omega}_{max}) \]

P2: prior não informativa em \(d\): \[P(d) \sim {\rm Unif}(0, +\infty)\]

ou, numa versão mais geral, fracamente informativo em \(d\): \[P(d) \sim {\rm Unif}(d_{min},d_{max})\]

Vamos examinar melhor o significado de um prior como \(P2\). Ele supõe que todas as distâncias são igualmente prováveis (num intervalo de 0 a \(\infty\) ou em um intervalo finito). Suponha que tenhamos um certo número de estrelas, \(N\), num ângulo sólido \(\Omega = A/d^2\). Conforme \(d\) aumenta, a área subtendida pelo ângulo sólido aumenta e a densidade de estrelas (\(N/A\)) diminui. Então, priores como esse, que supõem todas as distâncias igualmente prováveis, assumem implicitamente que existem menos estrelas por unidade de volume a grandes distâncias que a pequenas distâncias.

P3: prior informativo de Bailer-Jones et al. (2018) \[P(d) \propto d^2 e^{-d/L}\] para \(0 \le d_{min} < d < d_{max}\). Essa expressão corresponde a uma distribuição gama com fator de forma 3 e escala \(L\) (que é um hiperparâmetro).

Vamos ilustrar esses priores (unidades arbitrárias):

par(mfrow=c(1,1))
d = seq(0.1,10,0.1)
y1 = 1/d
y2 = 10*d^2*exp(-d/1)
# xaxt='n', yaxt='n': removem os valores nos gráficos
plot(d,y1,xlab='distância d', ylab='P(d)', type='l', col='salmon', xaxt='n', yaxt='n', lwd=3, main='prior')
abline(h=4,col='green',lwd=3)
lines(d,y2,lwd=3,col='blue')
legend('topright',legend=c('uniforme em paralaxe', 'uniforme em distância', 'Bailer-Jones et al. (2018)'), col=c('salmon','green','blue'),  text.col = c('salmon','green','blue'),bty='n')

3.4 Posterior

Queremos determinar a distância \(d\) a partir da medida da paralaxe \({\tilde \omega}\) de uma estrela e seu erro \(\sigma\) (dados \(D\)).

A função de distribuição de probabilidades da distância para esses dados é \[P(d|D) \propto P(D|d) P(d)\] e vai depender do prior: \[{\mathbf P1:} ~~~~ P(d|D) \propto {1 \over d} {\rm exp}\Big[ -{(\tilde{\omega} -1/d)^2 \over 2 \sigma^2} \Big]\] \[{\mathbf P2:} ~~~~ P(d|D) \propto {\rm exp}\Big[ -{(\tilde{\omega} -1/d)^2 \over 2 \sigma^2} \Big]\] \[{\mathbf P3:} ~~~~ P(d|D,L) \propto d^2 {\rm exp}\Big[ -{d \over L} -{(\tilde{\omega} -1/d)^2 \over 2 \sigma^2} \Big]\] para \(d_{min} < d < d_{max}\) (e 0 fora desse intervalo).

Vamos ver a forma desses posteriores (unidades arbitrárias):

par(mfrow=c(1,1))
d = seq(0.1,10,0.01)
omega = 1
sigma = 0.5
L=2
y1 = 1/d^2*exp(-(omega-1/d)^2/(2*sigma^2))
y2 = exp(-(omega-1/d)^2/(2*sigma^2))
y3 = d^2*exp(-d/L-(omega-1/d)^2/(2*sigma^2))
# xaxt='n', yaxt='n': removem os valores nos gráficos
plot(d,y1,xlab='distância d', ylab='P(d|D)', type='l', col='salmon', xaxt='n', yaxt='n', lwd=3, main='posterior')
lines(d,y2,lwd=3,col='green')
lines(d,y3,lwd=3,col='blue')
legend('topright',legend=c('uniforme em paralaxe', 'uniforme em distância', 'Bailer-Jones et al. (2018)'), col=c('salmon','green','blue'),  text.col = c('salmon','green','blue'),bty='n')

3.5 Exemplo prático

O aglomerado estelar das Hyades é um dos mais próximos de nós. O satélite Gaia mediu a distância de várias de suas estrelas e, para uma delas, obteve \(21.7321 \pm 0.0541\) mas.

# dados: a paralaxe e seu erro (em mas)
w = 21.7321
sigma = 0.0541

# solução de máxima verossimilhança (em pc):
dMV = 1/w*1000
dMV
## [1] 46.01488

Vamos calcular os posteriores para os 3 modelos acima:

# distâncias em pc
dpc = seq(45,47,0.01)
# e em kpc
d = dpc/1000

# posteriores (não-normalizados)
p1 = function(d){exp(-(w-1/d)^2/(2*sigma^2))/d}
p2 = function(d){exp(-(w-1/d)^2/(2*sigma^2))}
p3 = function(d){exp(-d/L-(w-1/d)^2/(2*sigma^2))*d^2}

# visualização
par(mfrow=c(3,1))
plot(dpc,p1(d)/max(p1(d)),xlab='distância d (pc)', ylab='P(d|D)', type='l', col='salmon', yaxt='n', lwd=3, main='posterior')
legend('topright',legend='uniforme em paralaxe', col='salmon',  text.col = 'salmon',bty='n')
abline(v=dMV)
plot(dpc,p2(d)/max(p2(d)),xlab='distância d (pc)', ylab='P(d|D)', type='l', col='green', yaxt='n', lwd=3)
legend('topright',legend='uniforme em distância', col='green',  text.col = 'green',bty='n')
abline(v=dMV)
plot(dpc,p3(d)/max(p3(d)),xlab='distância d (pc)', ylab='P(d|D)', type='l', col='blue', yaxt='n', lwd=3)
legend('topright',legend='Bailer-Jones', col='blue',  text.col = 'blue',bty='n')
abline(v=dMV)

Note que os posteriores são virtualmente idênticos! Isso acontece porque o erro na paralaxe é muito pequeno e a informação dos dados, não dos priores, domina.

Vamos determinar o intervalo de credibilidade do primeiro posterior. Para isso vamos, inicialmente, calcular sua função cumulativa. Pelas figuras acima podemos ver que o intervalo \(45 < d < 47\) pc contém o grosso das probabilidades.

d1 = 45/1000 # distâncias em kpc
d2 = 47/1000
F1 = function(d) integrate(p1,d1,d)$value/integrate(p1,d1,d2)$value

e agora vamos determinar o intervalo de credibilidade de 95%, isto é, a probabilidade da distância estar neste intervalo é 95%. Para isso vamos deteminar os valores de \(d\) que correspondem a 0.025 e 0.975 na função cumulativa:

# vamos fazer isso definindo uma função
q = 0.025
h = function(d){F1(d)-q}
# e obtendo d resolvendo a equação h=0 supondo que d está entre d1 e d2
dinf = uniroot(h, interval=c(d1,d2))$root
# solução em pc
round(dinf*1000,2)
## [1] 45.8
# idem para dsup
q = 0.975
# e obtendo d resolvendo a equação h=0 supondo que d está entre d1 e d2
dsup = uniroot(h, interval=c(d1,d2))$root
# solução em pc
round(dsup*1000,2)
## [1] 46.26

Assim, vemos que, pelo primeiro modelo, temos 95% de probabilidade de \(d\) estar entre 45.8 pc e 46.26 pc.

Vamos agora ver o que acontece se o erro fosse maior. A razão sinal-ruído dos dados atuais é \(\tilde{\omega}/\sigma\) ~ 400. E se fosse ~ 5?

Vamos supor que \(\tilde{\omega}\)=21.7 mas e \(\sigma=4.3\) mas e ver como ficam os posteriores:

# dados
wl = 21.7
sigmal = 4.3

# distâncias em pc
dpc = seq(0.1,100,0.01)
# e em kpc
d = dpc/1000

# posteriores
p1 = function(d){exp(-(wl-1/d)^2/(2*sigmal^2))/d}
p2 = function(d){exp(-(wl-1/d)^2/(2*sigmal^2))}
p3 = function(d){exp(-d/L-(wl-1/d)^2/(2*sigmal^2))*d^2}

# visualização
par(mfrow=c(3,1))
plot(dpc,p1(d)/max(p1(d)),xlab='distância d (pc)', ylab='P(d|D)', type='l', col='salmon', yaxt='n', lwd=3, main='posterior')
legend('topright',legend='uniforme em paralaxe', col='salmon',  text.col = 'salmon',bty='n')
abline(v=dMV)
plot(dpc,p2(d)/max(p2(d)),xlab='distância d (pc)', ylab='P(d|D)', type='l', col='green', yaxt='n', lwd=3)
legend('topright',legend='uniforme em distância', col='green',  text.col = 'green',bty='n')
abline(v=dMV)
plot(dpc,p3(d)/max(p3(d)),xlab='distância d (pc)', ylab='P(d|D)', type='l', col='blue', yaxt='n', lwd=3)
legend('topright',legend='Bailer-Jones', col='blue',  text.col = 'blue',bty='n')
abline(v=dMV)

Agora pode-se verificar que os diferentes priores levam a resultados diferentes.

Vamos comparar os intervalos de credibilidade desses posteriores.

# funções cumulativas- vamos considerar outro intervalo de integração
d1 = 20/1000 # distâncias em kpc
d2 = 200/1000
F1 = function(d) integrate(p1,d1,d)$value/integrate(p1,d1,d2)$value
F2 = function(d) integrate(p2,d1,d)$value/integrate(p2,d1,d2)$value
F3 = function(d) integrate(p3,d1,d)$value/integrate(p3,d1,d2)$value

e agora vamos determinar o intervalo de credibilidade de 95%:

# para P1
q = 0.025
h = function(d){F1(d)-q}
# e obtendo d resolvendo a equação h=0 supondo que d está entre d1 e d2
dinf = uniroot(h, interval=c(d1,d2))$root
# solução em pc
dinf1=round(dinf*1000,2)

# para P2
h = function(d){F2(d)-q}
dinf = uniroot(h, interval=c(d1,d2))$root
dinf2=round(dinf*1000,2)

# para P3
h = function(d){F3(d)-q}
dinf = uniroot(h, interval=c(d1,d2))$root
dinf3=round(dinf*1000,2)

# idem para dsup
q = 0.975
h = function(d){F1(d)-q}
dsup = uniroot(h, interval=c(d1,d2))$root
dsup1 = round(dsup*1000,2)

h = function(d){F2(d)-q}
dsup = uniroot(h, interval=c(d1,d2))$root
dsup2 = round(dsup*1000,2)

h = function(d){F3(d)-q}
dsup = uniroot(h, interval=c(d1,d2))$root
dsup3 = round(dsup*1000,2)

# resultados:
c('P1: ',dinf1,dsup1)
## [1] "P1: "  "34.06" "82.86"
c('P2: ',dinf2,dsup2)
## [1] "P2: "  "35.03" "93.96"
c('P3: ',dinf3,dsup3)
## [1] "P3: "   "37.55"  "139.69"

Vemos que, nesse caso, os intervalos de credibilidade dos 3 modelos fica bem diferente.

3.6 A evidência

Vamos usar a evidência para comparar os 3 modelos acima.

Para a estimativa dos parâmetros na seção anterior os priores não precisavam estar normalizados, mas para o cálculo da evidência isso é necessário.

Vamos supor que, a priori, supomos que a distância à estrela esteja entre \(d_{min} = 10\) pc e \(d_{max} = 100\) pc (ou a paralaxe esteja entre \(\tilde{\omega}_{min} = 1/d_{max} = 10\) mas e \(\tilde{\omega}_{max} = 1/d_{min} = 100\) mas).

Nesse caso temos:

dmin = 10 #pc
dmax = 100 #pc

aux = function(d) 1/d
cte = integrate(aux,dmin,dmax)$value
prior1 = function(d) 1/(d*cte)
prior2 = function(d) 1/(dmax-dmin)
prior3 = function(d){p3(d)/integrate(p3,dmin,dmax)$value}

e a evidência desses modelos é, então,

f1 = function(d) p1(d)*prior1(d)
e1 = integrate(f1,dmin,dmax)$value

f2 = function(d) p2(d)*prior2(d)
e2 = integrate(f2,dmin,dmax)$value

f3 = function(d) p3(d)*prior3(d)
e3 = integrate(f3,dmin,dmax)$value

c(e1,e2,e3)
## [1] 1.230461e-07 3.040767e-06 1.161636e-06

O modelo com maior evidência é o 2. Vamos comparar os modelos 1 e 3 com 2 (logaritmo natural dos odds):

log(e2/e1)
## [1] 3.207306
log(e2/e3)
## [1] 0.9622804

A confiança na diferença entre os modelos 2 e 1, de acordo com a escala de Jeffreys (modificada por Trotta), é moderada, mas entre 2 e 3 é fraca. O modelo 2 parece ser um pouco melhor que o 1, mas sua performance é semelhante ao 3.

Exercícios

  1. (Trotta) 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, reportou \(1.98 \pm 0.16\) seg. de arco.

Use a comparação bayesiana de hipóteses para saber se os dados favorecem Newton ou Einstein:

  • qual é o prior de Newton e o de Einstein para \(\alpha\)?
  • se as observações podem ser descritas por \(N(\mu,\sigma)\), escreva a verossimilhança dos dados e o fator de Bayes entre Einstein e Newton
  • compare os fatores de Bayes obtidos com as observações de Eddington e as de Crommelin. Qual expedição deu resultados de maior evidência?
  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.