aula de hoje:

  1. inferência bayesiana: análise sequencial ou conjunta?
  1. quando os métodos bayesianos e frequentistas divergem?
  1. o modelo beta-binomial
  1. o problema do farol
  1. um problema de astrofísica de altas energias
  1. quantos peixes tem no lago?
  1. estimativa dos parâmetros de uma gaussiana por ABC
  1. a estratégia de marketing da Swedish Fish Incorporated

Exercícios

1. análise de dados bayesiana: sequencial ou conjunta?

  • temos um conjunto de dados \(D=\{x_i\}\) que queremos modelar com uma distribuição normal, \(N(\mu,\sigma)\), onde conhecemos \(\sigma\) mas não conhecemos \(\mu\)
  • vamos usar o teorema de Bayes para estimar \(\mu\) dos dados. Para isso precisamos definir a verossimilhança dos dados e o prior dos parâmetros
  • verossimilhança dos dados: \[ P(D|\mu,\sigma) = \prod_{i=1}^n P(x_i|\mu,\sigma) \propto \exp \Bigg[ -\frac{1}{2\sigma^2}\sum_{i=1}^n(x_i-\mu)^2 \Bigg] \]
  • vamos adotar um prior também gaussiano para \(\mu\): \(P(\mu) \sim N(\mu_0,\tau_0)\)
  • nesse caso o posterior também é gaussiano: \(N( \mu_n, \tau_n)\) \[ \mu_n = \frac{\frac{\mu_0}{\tau_0^2}+\frac{n \bar x}{\sigma^2}}{\frac{1}{\tau_0^2}+\frac{n}{\sigma^2}}~~~~~~~~~~~ \frac{1}{\tau_n^2} = \frac{1}{\tau_0^2}+\frac{n}{\sigma^2}~~~~~~~~~~~ \bar x = \frac{1}{n}\sum_i^n x_i \]
  • simulação: vamos supor dados distribuídos de acordo com o modelo \(N(\mu=3,\sigma=1)\) e adotar o prior \(N( \mu_0 = -1,\tau_0=1)\), propositalmente bem viesado!

Vamos inicialmente visualizar o modelo e o prior:

# sequência de pontos entre -4 e 6 separados por 0.1
x = seq(-4, 6, by = .1) 

# vetor com os valores do modelo gaussiano para cada valor de x
y = dnorm(x,mean=3,sd=1) 

# visualização do prior e do modelo
# modelo
plot(x,y,main="prior (V) e modelo (A)",type="l",col='blue',lwd=3)  

# prior
lines(x,dnorm(x,mean=-1,sd=1),col='red',lwd=3)  # prior

Vamos gerar 1000 dados com o modelo \(N(\mu=3,\sigma=1)\) e ver como o posterior muda com o número de dados.

# semente do gerador de números aleatórios para assegurar a reprodutibilidade dos resultados
set.seed(123)

# parâmetros do modelo
n=1000
musim=3
sigma=1

# parâmetros do prior
mu0=-1
tau0=1

# simulação dos dados
d=rnorm(n,mean=musim,sd=sigma)

# janela 3 x 2 para as figuras:
par(mfrow=c(3,2))

#gaussianas normalizadas pelo máximo para facilitar visualização
plot(x,y/max(y),main="prior (vermelho) e modelo (azul)",type="l",col='blue',lwd=3)  
yy=dnorm(x,mean=mu0,sd=tau0)
lines(x,yy/max(yy),col='red',lwd=3)

# loop para fazer os gráficos para valores escolhidos do número de pontos, m:
for (m in c(1,2,10,100,1000)){
  
# parametros do posterior em cada caso:
xmedio = mean(d[1:m])  
taun = 1/sqrt(1/tau0^2 + m/sigma^2)
mun = (mu0/tau0^2+m*xmedio/sigma^2)/(1/tau0^2 + m/sigma^2)

#valores do posterior
yyy= dnorm(x,mean=mun,sd=taun)

# plot do posterior
plot(x,y/max(y),type="l",col='blue', main=sprintf(paste('m =  ',m,';  posterior em verde')),lwd=3)
lines(x,yy/max(yy),col='red',lwd=3)
lines(x,yyy/max(yyy),col='green',lwd=3)
}

Note que:

  • a moda do posterior migrou rapidamente de próximo da moda do prior para a dos dados simulados
  • a largura do posterior diminui conforme o número de dados aumenta

2. 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)=(1-p)^3\)
  • vamos supor o prior \(P(p)\) uniforme entre 0 e 1
  • verossimilhança \(P(D|p) \propto p^5(1-p)^3\)

logo, \[ 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:

#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

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

3. o modelo beta-binomial

A distribuição binomial descreve problemas onde os dados são do tipo 0/1, TRUE/FALSE, Detectado/não-Detectado, Cara/Coroa; em geral sucesso/não-sucesso

A distribuição binomial tem um único parâmetro, \(\theta\): a probabilidade de sucesso em cada tentativa

Probabilidade de \(n\) sucessos em \(N\) tentativas independentes: \[ P(n|\theta,N)=\binom Nn \theta^n (1-\theta)^{N-n} = {N! \over n!(N-n)!} \theta^n (1-\theta)^{N-n}\] média e variância: \[ \bar n = \sum_{n=0}^N n P(n) = N\theta ~~~~~~ \sigma^2 = \sum_{n=0}^N (n-\bar n)^2 P(n) = N\theta(1-\theta) \]

Dado um conjunto de dados, vamos fazer um modelo bayesiano para estimar \(\theta\) (\(0 \le \theta \le 1\))

Vamos escrever a verossimilhança como a distribuição binomial \[ P(\theta|n,N) = \binom{N}{n} \theta^n(1-\theta)^{N-n} \]

e o prior como uma função Beta, já que a variável \(\theta\) é definida entre 0 e 1: \[ Beta(\theta|a,b) = \Bigg[ \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}\Bigg] \theta^{a-1} (1-\theta)^{b-1} \] \[(a,b ~ > 0) \]

posterior: \[ P(\theta|D) \propto \theta^{n+a-1}(1-\theta)^{N-n+b-1} \] \[ \propto Beta(\theta|n+a,N-n+b) \]

note que, nesse caso, tanto o prior quanto o posterior são funções Beta: Beta é um prior conjugado da distribuição binomial

Exemplo: baseado em https://www.r-exercises.com/2017/12/03/basic-bayesian-inference-for-mcmc-techniques-exercises-part-1/

Um modelo de formação de galáxias prevê que em aglomerados de galáxias ricos a maior parte das galáxias luminosas sejam elípticas ou lenticulares (em baixos redshifts). Queremos analisar se esse é realmente o caso em uma amostra de aglomerados. Este problema é uma adaptação do descrito no link acima.

Seja \(\theta\) a fração de galáxias early type (ET: E e S0). O estudo de uma grande simulação revela que em uma amostra de \(N=1000\) galáxias luminosas em aglomerados ricos, \(n=735\) são desse tipo. Vamos classificar as galáxias como \(y=1\) se forem ET e \(y=0\) se não o forem.

Vamos ver, inicialmente, qual é a distribuição de probabilidade dos dados, plotando a verossimilhança:

# reprodutibilidade
set.seed(123)

plot(dbinom(x = seq(1, 1000, 1), size = 1000, prob = 735/1000),
     type = "l", lwd = 3, main = "verossimilhança", ylab = "densidade de probabilidades", xlab = "número de galáxias early type", col='blue')
abline(v = 735, col = "red",lwd=2)

Vamos considerar inicialmente um prior uniforme para \(\theta\); ele corresponde a \(Beta(\theta|a=1,b=1)\)

a=1
b=1
x=seq(0,1,0.001)
y = dbeta(x,a,b)
plot(x,y,type='l',col='blue', main='prior', ylab = "densidade de probabilidades",xlab=expression(theta),lwd=3)

Logo, o posterior será \[ P(\theta|D,a,b) \propto Beta(\theta|a=736,b=266):\]

a=1
b=1
n=735
N=1000
y = dbeta(x,n+a,N-n+b)
plot(x,y,type='l',col='green', main='posterior', ylab = "densidade de probabilidades",xlab=expression(theta),lwd=3)

# moda   
# https://stackoverflow.com/questions/2547402/how-to-find-the-statistical-mode
Mode <- function(x) {
  ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))]
}

M=100000
post_sample = rbeta(M, shape1 = n+a, shape2 = N-n+b)
Mode(post_sample)
## [1] 0.7433852
# intervalo de de credibilidade de 95%:
quantile(post_sample, probs = c(0.025, 0.975))
##      2.5%     97.5% 
## 0.7067023 0.7614554

Uma outro estudo indica que apenas 60+/-20% das galáxias luminosas são ET. Queremos usar esta informação como prior. Como podemos fazer isso?

A distribuição beta tem dois hiperparâmetros ("shape parameters") que podem ser obtidos da média e variância da distribuição:

estBetaParams <- function(mean, variance) {
  alpha <- ((1 - mean) / variance - 1 / mean) * mean^2
  beta <- alpha * (1 / mean - 1)
  return(params = list(a = alpha, b = beta))
}

Assim, em nosso exemplo o prior de \(\theta\) tem parâmetros:

prior.params <- estBetaParams(mean = 0.6, variance = 0.2^2)
prior.params
## $a
## [1] 3
## 
## $b
## [1] 2

Vamos representar graficamente este prior:

x_grid <- seq(from = 0,to = 1,  by = 0.001)
# note que seq deve estar no intervalo de Beta, [0,1]

plot(x = x_grid, dbeta(x = x_grid, shape1 = prior.params$a, shape2 = prior.params$b), lwd = 3, col = 'blue', ylab = "densidade", xlab = expression(theta), type = "l", main = "prior")

e o posterior será:

N <- 1000
n <- 735
# hiper-parâmetros do posterior
alpha_post <- prior.params$a + n
beta_post <- prior.params$b + N - n
#plot do posterior
plot(x = x_grid, dbeta(x = x_grid, shape1 = alpha_post, shape2 = beta_post), lwd = 3, type = "l", ylab = "densidade", xlab = expression(theta), col = "red", main = "posterior e prior")  
# prior
lines(x = x_grid, dbeta(x = x_grid, shape1 = prior.params$a, shape2 = prior.params$b), lwd = 3) 
legend("topright",c("Posterior","Prior"),col=c("red", "black"),lwd=c(3,3))
# valor de teta da verossimilhança
abline(v = 0.735, col = "green")

# moda
M=100000
post_sample = rbeta(M, shape1 = n+a, shape2 = N-n+b)
Mode(post_sample)
## [1] 0.7396033
# intervalo de credibilidade de 95%:
post_sample = rbeta(M, shape1 = prior.params$a, shape2 = prior.params$b)
quantile(post_sample, probs = c(0.025, 0.975))
##      2.5%     97.5% 
## 0.1942310 0.9326187

Pode-se notar que, nesse exemplo, a forma do prior afeta muito pouco o posterior. Compare a moda e o intervalo de confiança de \(\theta\) nos dois casos. Se o número de dados fosse bem menor as diferenças seriam maiores!

4.o problema do farol

Este problema é atribuído a Gull (1988) e discutido em Data Analysis - A Bayesian Tutorial, Sivia & Silling (2a. ed., 2006).

Um farol está em uma ilha em uma posição \(\alpha\) ao longo de uma costa reta e a uma distância \(\beta\) no mar; girando, ele emite uma série de pulsos curtos altamente colimados em intervalos de tempo aleatórios (e portanto em azimutes também aleatórios); \(N\) destes pulsos são detectados por sensores na costa, mas só as posições \(\{ x_k \}\), não as direções.

Onde está o farol? Quanto valem \(\alpha\) e \(\beta\)?

Dado um conjunto de medidas, \(x\), o posterior de \(\alpha\) e \(\beta\) é \[ P(\alpha,\beta|x) \propto P(x|\alpha,\beta) P(\alpha,\beta)\]

Vimos que os dados devem seguir uma distribuição de Cauchy, que é a verossimilhança do problema: \[P(x|\alpha,\beta) = {\beta \over \pi [\beta^2 + (x-\alpha)^2]}. \]

Vamos supor priores não informativos para \(\alpha\) e \(\beta\): \(P(\alpha,\beta) =\) const.

Nesse caso o posterior é proporcional à verossimilhança.

Vamos simular dados para este problema. Vamos supor que \(\theta\) está uniformemente distribuído entre \(-\pi/2\) e \(\pi/2\) e vamos simular valores de \(x\) via \[ x_k = \alpha_0 + \beta_0 \tan \theta_k, \] onde \((\alpha_0,\beta_0)\) são os valores verdadeiros dos parâmetros do modelo.

simulação dos dados:

# reprodutibilidade
set.seed(123)

#parâmetros do modelo
npar = 2
alpha0=0.5
beta0=1.5

# número de observações
npts=100

# simulação das observações:
# simulação de theta, em radianos
theta=runif(npts,min=-pi/2,max=pi/2)

# os 10 primeiros valores de theta simulados, em graus:
theta[1:10]*180/pi
##  [1] -38.236046  51.894924 -16.384154  68.943133  79.284111 -81.799830
##  [7]   5.058988  70.635428   9.258303  -7.809348
# simulação de x
x=alpha0+beta0*tan(theta)

# distribuição das observações:
par(mfrow = c(1,2))
# todas as observações
hist(x,col='salmon',breaks=1000,ylab='frequência',main='')
# observações próximas ao máximo
hist(x,col='salmon',breaks=1000,xlim=c(-10,10),ylab='frequência',main='')

Note que a distribuição observada de \(x\) contém valores muito grandes (em módulo): os outliers de \(x\) são típicos de uma distribuição de Cauchy. Vejam os valores extremos nesse caso:

summary(x)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -763.7199   -1.0433    0.3409   -6.1294    2.0526   83.8149

Comparem o desvio padrão e o desvio absoluto mediano:

sd(x)
## [1] 77.14604
mad(x)
## [1] 2.249799

notem que o desvio padrão é mais sensível a outliers que o MAD.

Para determinar os parâmetros vamos definir nosso modelo probabilístico.

Vamos começar definindo uma função para determinar o log da verossimilhança de uma distribuição de Cauchy, de um único dado. Em muitos casos a verossimilhança dos dados pode ter valores muito pequenos (por exemplo \(10^{-100}\)), e quando muitos valores assim são multiplicados, podem aparecer problemas numéricos (como "underflow"), que abortam a execução de um programa ou comprometem a precisão dos resultados. O uso de logaritmos evita isso.

verossimilhança dos dados (em log):

# como as medidas são independentes, o log da verossimilhança é a soma dos logs das verossimilhanças de cada dado individualmente

# log da verossimilhança dos dados
# o input é o vetor de parâmetros: param[1] = alfa e param[2] = beta    
logver = function(param) sum(log(param[2]/pi/(param[2]^2+(x-param[1])^2) ))

prior (em log):

# log do prior
logprior = function(param){ 
  logpriora = dunif(param[1], min=-50, max=50, log = T)
  logpriorb = dunif(param[2], min=0, max=50, log = T)
  return(logpriora+logpriorb)}

posterior (em log):

logposterior = function(param){logver(param) + logprior(param)}

função proposta: \(N(0,0.1)\) para cada parâmetro

funcao_proposta = function(param){
    rnorm(npar,mean = param, sd= rep(0.1,npar))}

Vamos usar a função mcmc, aqui ajustada para um posterior em log. A função recebe como valores de entrada um chute inicial para os valores dos parâmetros e o número de iterações do algoritmo, e retorna uma cadeia de valores de parâmetros:

mcmc = function(xini, niter){
    cadeia = array(dim = c(niter+1,npar))
    cadeia[1,] = xini
    for (i in 1:niter){
        proposta = funcao_proposta(cadeia[i,])
        probab = exp(logposterior(proposta) - logposterior(cadeia[i,]))
        ifelse ((runif(1) < probab),
            yes = (cadeia[i+1,] = proposta),
            no = (cadeia[i+1,] = cadeia[i,]))
    }    
    return(cadeia)
}

simulacão de MCMC:

# simulacao
xinicial = c(2,2)
nsim=10000
simulacao = mcmc(xinicial, nsim)

# vamos supor um burn-in de 1000 iterações:
# (o objetivo do burn-in é remover as simulações iniciais, que podem estar fora da região de maior probabilidade dos parâmetros)
burnin=1000

# taxa de aceitação
aceitacao = 1-mean(duplicated(simulacao[-(1:burnin),]))
aceitacao
## [1] 0.7596934

visualização da cadeia:

# visualização, sem o burn-in
# a linha vermelha é o valor médio do parâmetro
par(mfrow = c(2,2))
hist(simulacao[-(1:burnin),1],breaks=30, main=expression(paste("Posterior de ",alpha)), xlab=expression(alpha),ylab='frequência',col='gold')
abline(v = mean(simulacao[-(1:burnin),1]),lwd=2)
hist(simulacao[-(1:burnin),2],breaks=30, main=expression(paste("Posterior de ",beta)), xlab=expression(beta),ylab='frequência',col='gold')
abline(v = mean(simulacao[-(1:burnin),2]),lwd=2)
plot(simulacao[-(1:burnin),1], type = "l", xlab="iter" , main = expression(paste("cadeia de ",alpha)),ylab='frequência',col='gray50')
abline(h = mean(simulacao[-(1:burnin),1]), col="red",lwd=2)
plot(simulacao[-(1:burnin),2], type = "l", xlab="iter" , main =  expression(paste("cadeia de ",beta)),ylab='frequência',col='gray50')
abline(h = mean(simulacao[-(1:burnin),2]), col="red",lwd=2)

estatísticas dos parâmetros: compare os resultados com os parâmetros usados na simulação, \((\alpha_0,\beta_0) = (0.5,1.5)\)

# valores médios:
ca=simulacao[-(1:burnin),1]
cb=simulacao[-(1:burnin),2]
media.a=mean(ca)
media.a
## [1] 0.4386645
media.b=mean(cb)
media.b
## [1] 1.516981
# desvios padrão
sd(ca)
## [1] 0.2203618
sd(cb)
## [1] 0.2171015
# intervalo de credibilidade dos parâmetros:
quantile(ca, c(0.025, 0.975))
##      2.5%     97.5% 
## 0.0167786 0.8812929
quantile(cb, c(0.025, 0.975))
##     2.5%    97.5% 
## 1.130389 1.967030

visualização do posterior:

library(KernSmooth)
## KernSmooth 2.23 loaded
## Copyright M. P. Wand 1997-2009
# mapa de densidades
par(mfrow = c(1,1))
smoothScatter(ca,cb,nrpoints=0,add=FALSE,xlab=expression(alpha),ylab=expression(beta))
points(media.a,media.b,pch=3,col="red")   
points(alpha0,beta0,pch=3,col="green")  

# contornos
Ndim = 101
z = cbind(ca,cb)
xmax = max(ca)
xmin = min(ca)
ymax = max(cb)
ymin = min(cb)
deltax = (xmax-xmin)
deltay = (ymax-ymin)
bw.x = 6*deltax/Ndim
bw.y = 6*deltay/Ndim
dens = bkde2D(z,bandwidth=c(bw.x,bw.y),gridsize=c(Ndim,Ndim)) 
contour(dens$x1,dens$x2,dens$fhat,add=T,drawlabels=F,nlevels=5)

# vermelho: parâmetros médios da simulação       verde: parâmetros verdadeiros

5. um problema de astrofísica de altas energias

Um telescópio de raios-X coleta fótons em diversos intervalos (bins) de energia. A tabela a seguir contém a energia média e o número de contagens em cada bin:

#dados
dados = read.table(file="/home/laerte/cursos/dados/2022/dados/contagensX.dat",header=TRUE)
E = dados[,1]
N = dados[,2]

#visualização:
plot(E,N,type='h',lwd=2,col='blue')

Vamos supor que as contagens no bin \(i\) (correspondente à energia \(E_i\)) podem ser modeladas como um processo poissoniano, \(N_i \sim Poisson(C_i)\), onde \(C_i\) é o número esperado de fótons no bin, que vamos modelar como \[C_i = \alpha E_i^{-\beta},\] onde \(\alpha\) e \(\beta\) são os parâmetros do modelo.

Vamos fazer um modelo bayesiano para este processo, assumindo priores uniformes:

#reprodutibilidade
set.seed(123)

# log da verossimilhança dos dados
param=NULL
logver = function(param){
NC=param[1]*E^(-param[2])
return(-0.5*sum((N-NC)^2))
}

# log do prior
logprior = function(param){ 
  logpriora = dunif(param[1], min=20, max=300, log = T)
  logpriorb = dunif(param[2], min=1.5, max=3, log = T)
return(logpriora+logpriorb)}

#posterior (em log):
logposterior = function(param){logver(param) + logprior(param)}

#função proposta: 
funcao_proposta = function(param){rnorm(npar,mean = param, sd= c(1,0.1))}

#mcmc
mcmc = function(xini, niter){
    cadeia = array(dim = c(niter+1,npar))
    cadeia[1,] = xini
    for (i in 1:niter){
        proposta = funcao_proposta(cadeia[i,])
        probab = exp(logposterior(proposta) - logposterior(cadeia[i,]))
        ifelse ((runif(1) < probab),
            yes = (cadeia[i+1,] = proposta),
            no = (cadeia[i+1,] = cadeia[i,]))
    }    
    return(cadeia)
}

#simulacão de MCMC:
# simulacao
xinicial = c(100,3)
nsim=1000000
simulacao = mcmc(xinicial, nsim)

# vamos supor um burn-in de 1000 iterações:
burnin=1000

# taxa de aceitação
aceitacao = 1-mean(duplicated(simulacao[-(1:burnin),]))
aceitacao
## [1] 0.1640789
#sumário dos resultados, sem o burn-in
p1 = simulacao[-(1:burnin),1]
p2 = simulacao[-(1:burnin),2]
summary(p1)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   170.3   200.9   207.1   207.4   213.7   246.1
summary(p2)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.366   2.540   2.573   2.573   2.606   2.785
# intervalo de de credibilidade de 95%:
quantile(p1, probs = c(0.025, 0.975))
##     2.5%    97.5% 
## 188.9543 227.6866
quantile(p2, probs = c(0.025, 0.975))
##     2.5%    97.5% 
## 2.475154 2.673941
# visualização, sem o burn-in
# a linha vermelha é o valor médio do parâmetro
par(mfrow = c(2,2))
hist(simulacao[-(1:burnin),1],breaks=30, main=expression(paste("Posterior de ",alpha)), xlab=expression(alpha),ylab='frequência',col='gold')
abline(v = mean(simulacao[-(1:burnin),1]),lwd=2)
hist(simulacao[-(1:burnin),2],breaks=30, main=expression(paste("Posterior de ",beta)), xlab=expression(beta),ylab='frequência',col='gold')
abline(v = mean(simulacao[-(1:burnin),2]),lwd=2)
plot(simulacao[-(1:burnin),1], type = "l", xlab="iter" , main = expression(paste("cadeia de ",alpha)),ylab='frequência',col='gray50')
abline(h = mean(simulacao[-(1:burnin),1]), col="red",lwd=2)
plot(simulacao[-(1:burnin),2], type = "l", xlab="iter" , main =  expression(paste("cadeia de ",beta)),ylab='frequência',col='gray50')
abline(h = mean(simulacao[-(1:burnin),2]), col="red",lwd=2)

# mapa de densidades
par(mfrow = c(1,1))
smoothScatter(simulacao[-(1:burnin),1],simulacao[-(1:burnin),2],nrpoints=0,add=FALSE,xlab=expression(alpha),ylab=expression(beta))
points(mean(p1),mean(p2),pch=3,col="red",lwd=3)   

# ajuste com os parâmetros médios:
par(mfrow=c(1,1))
NC=mean(p1)*E^(-mean(p2))
plot(E,N,type='h',lwd=2,col='blue')
lines(E,NC,lwd=3,col='salmon')

6. quantos peixes tem no lago?

Vamos a um lago e pescamos 7 peixes. Marcamos cada um e os devolvemos ao lago. Alguns dias depois voltamos ao lago, pescamos 20 peixes (dia de sorte!) e verificamos que 3 deles estão marcados. Quantos peixes tem no lago?

Seja \(T\) o número total de peixes no lago, \(M\) o número de peixes marcados, \(K\) o número de peixes pescados e \(X\) o número dentre eles que estão marcados.

Dados: \(M=7\), \(K=20\) e \(X=3\); queremos determinar \(T\)

Solução de MV: se supomos que a fração de peixes marcados (x/K) resultante de nossa medida (pescaria) é representativa, então \[ \hat T = MK/X = 7*20/3 \simeq 46.67\]

Pode-se mostrar que a verossimilhança é dada pela distribuição multinomial: \[ P(X) = \frac{{M\choose X}{T-M\choose K-X}}{{T\choose K}} \]

Para determinar o posterior de \(T\), vamos supor um prior uniforme entre 20 (o número mínimo de peixes no lago) e 400

# reprodutibilidade
set.seed(123)

# dados
X=3
M=7
K=20
# solução de máxima verossimilhança
MV = M*K/X
MV
## [1] 46.66667
# como o prior é uniforme, o posterior de T é igual a verossimilança
# vamos escrevê-lo como uma função:
prob = function(T) {choose(M,X)*choose(T-M,K-X)/choose(T,K)}
#visualização do posterior
T=seq(20,400,1)
y=prob(T)
plot(T,y,col='salmon',pch=20,type='h',lwd=2,xlim=c(0,410))
# valor de MV
abline(v=MV,lwd=2)

# intervalo de credibilidade de 95%
s= cumsum(y)/sum(y)
T025 = 20+max(which(s <= 0.025))
T975 = 20+min(which(s >= 0.975))
c(T025,T975)
## [1]  32 289
abline(v=T025,col='blue')
abline(v=T975,col='blue')

Pode-se notar como a cauda da distribuição é extensa!

Vamos agora resolver o mesmo problema simulando modelos gerativos e usando ABC

adaptado de https://rpubs.com/rasmusab/live_coding_user_2015_bayes_tutorial

# dados:
# número de peixes marcados:
M = 7
# número de peixes pescados: 
K = 20
# número de peixes marcados entre os pescados: 
X = 3

# número de simulações
nsim <- 100000

# vamos considerar um prior uniforme entre 20 e 400 e amostrar o número de peixes com esse prior
# n_peixes: número de peixes no lago
n_peixes <- sample(20:400, nsim, replace = TRUE)

# aqui definimos o modelo gerativo:
# para cada um dos n_peixes, 0: não marcado, 1: marcado
# crio uma sequência de 0s seguida de M pontos marcados

# pescaria: número de peixes pescados:
pesca_peixes <- function(n_peixes) {
peixes <- rep(0:1, c(n_peixes - M, M))
sum(sample(peixes, K))
}

# número de peixes marcados em cada simulação
n_marcados <- rep(NA, nsim)
for(i in 1:nsim) {n_marcados[i] <- pesca_peixes(n_peixes[i])}

# posterior por ABC:
#seleciona-se as simulações onde se obtém o número de peixes marcados correto 
post_npeixes <- n_peixes[n_marcados == X]

hist(post_npeixes,breaks=100,col='salmon',main='posterior',xlab='no. de peixes')
abline(v=MV,lwd=3)

# intervalo de credibilidade de 95%
quantile(post_npeixes, c(0.025, 0.975))
##  2.5% 97.5% 
##    31   298
abline(v=T025,col='blue')
abline(v=T975,col='blue')

# compare com o intervalo obtido acima

7. Estimativa dos parâmetros de uma gaussiana por ABC

Vamos amostrar N=10 valores de uma gaussiana de média 2 e desvio padrão 1 e usar ABC para inferir os parâmetros da gaussiana.

Vamos criar um modelo gerativo e usar a média e a variância dos dados como as estatísticas relevantes.

Vamos usar também um valor máximo eps para a diferença dessas estatisticas entre os dados e os modelos gerativos para fazer a aceitação por ABC.

Vamos considerar que temos 2 parâmetros, 1) a média e 2) o desvio padrão.

# reprodutibilidade
set.seed(123)

# parâmetros da simulação: gaussiana de média 2 e desvio padrão 1
npar = 2
parsim = rep(0,npar)
parsim[1] = 2   # média
parsim[2] = 1   # desvio padrão

# dados simulados:
N = 10
dados = rnorm(N, mean = parsim[1], sd = parsim[2])

# estatísticas dos dados:
pard = NULL
pard[1] = mean(dados)
pard[2] = sd(dados)

# para uma amostra de dados ser aceita por ABC, sua média e desvio padrão devem diferir da média e desvio padrão dos dados por no máximo eps:
eps_media = 0.1
eps_sd =  0.2

# função que examina se um conjunto de dados satisfaz (T) ou não (F) as condições de ABC
aceitacao_ABC = function(par){
# o desvio padrão tem que ser positivo
  if (par[2] <= 0) return(F)
  
# modelo gerativo dos dados:
amostras = rnorm(N, mean =par[1], sd = par[2])
delta_media = abs(mean(amostras)-pard[1])
delta_sd = abs(sd(amostras)-pard[2])
ifelse(delta_media < eps_media & delta_sd < eps_sd, yes=T, no= F)
}

# vamos agora usar MCMC usando a aceitação para ABC:
MCMC_ABC <- function(pini, iteracoes){
    cadeia = array(dim = c(iteracoes+1,npar))
    cadeia[1,] = pini
    for (i in 1:iteracoes){
        # função proposta:
        proposta = rnorm(npar,mean = cadeia[i,], sd= c(0.5,0.5))
        # verifica se a proposta é aceita ou não
        if(aceitacao_ABC(proposta)){    
            cadeia[i+1,] = proposta
        }else{
            cadeia[i+1,] = cadeia[i,]
        }
    }
    return(cadeia)
}

niter = 500000
parini = c(1.5,1.5)

posterior <- MCMC_ABC(parini,niter)
medias = posterior[,1]
desvp = posterior[,2]

# como estamos fazendo muitas simulações o burn-in não é importante

par(mfrow = c(2, 2))
it = 1:nrow(posterior)
plot(it,medias,type='l',col='gray',main='',xlab='iteração',ylab='média')
abline(h=mean(medias),lwd=2)
abline(h=parsim[1],lwd=2,col='red')
hist(medias,col='lightblue1',main='')
abline(v=mean(medias),lwd=2)
abline(v=parsim[1],lwd=2,col='red')
plot(it,desvp,type='l',col='gray',main='',xlab='iteração',ylab='desvio padrão')
abline(h=mean(desvp),lwd=2)
abline(h=parsim[2],lwd=2,col='red')
hist(desvp,col='lightblue1',main='')
abline(v=mean(desvp),lwd=2)
abline(v=parsim[2],lwd=2,col='red')

# valores estimados e intervalos de credibilidade de 95%
c(mean(medias),quantile(medias, c(0.025, 0.975)))
##              2.5%    97.5% 
## 2.062755 1.325445 2.764124
c(mean(desvp),quantile(desvp, c(0.025, 0.975)))
##                2.5%     97.5% 
## 1.1126283 0.6318008 1.9837909

Compare os valores obtidos por ABC com os iniciais, usados para simular os dados.

8. Teste bayesiano para a estratégia de marketing da Swedish Fish Incorporated

Este exercício é uma adaptação do blog de Rasmus Baath: https://www.sumsar.net/files/posts/2017-bayesian-tutorial-exercises/modeling_exercise1.html

A Swedish Fish Incorporated (SFI) é a maior empresa sueca de entrega de peixe por correspondência. Eles agora estão tentando entrar no lucrativo mercado dinamarquês vendendo assinaturas de um ano de salmão.

O Departamento de Marketing fez um estudo piloto e experimentou uma estratégia de marketing, chamada de método A, enviando um e-mail com um folheto colorido que convida as pessoas a se inscreverem para uma assinatura de salmão por um ano.

O Departamento de Marketing enviou 16 correspondências do tipo A. Seis dinamarqueses que receberam a correspondência assinaram por um ano de salmão e o Departamento de Marketing agora quer saber: quão bom é o método A?

No final dessa seção você encontrará uma solução. Mas tente você mesmo primeiro!

Questão I) Construa um modelo bayesiano que responda à pergunta: qual seria a taxa de adesão se o método A fosse usado em um número maior de pessoas?

Dica 1: A resposta não é um único número, mas uma distribuição sobre as taxas de adesão prováveis.

Dica 2: Como parte de seu modelo gerativo, você desejará usar a distribuição binomial, que você pode amostrar em R usando rbinom(n, size = m, prob). A distribuição binomial simula o seguinte processo n vezes: o número de “sucessos” em m tentativas, onde a probabilidade de “sucesso” em cada tentativa é prob.

Dica 3: Um prior comumente usado para a probabilidade desconhecida de sucesso em uma distribuição binomial é uma distribuição uniforme de 0 a 1. Você pode amostrar esta distribuição executando runif(1, min = 0, max = 1)

Dica 4: Aqui está um esquema de código que você pode melhorar:

## reprodutibilidade
# set.seed(numero_USP)

## número de amostras extraídas do prior de prob
# n_amostras <- 10000

## aqui você obtém n_amostras amostras do prior
# prior <- ...  

## é sempre bom olhar o prior para ver se ele parece ok
# hist(prior) 

## aqui você define seu modelo gerativo, onde você obtém n_amostras em função do valor de prob
# modelo_gerativo <- function(parametros) {
# ...
#}

## aqui você simula dados usando os parâmetros do prior e o modelo gerativo 
# sim_dados <- rep(NA, n_amostras)
#for(i in 1:n_amostras) {
#  sim_dados[i] <- modelo_gerativo(prior[i])
#}

## aqui você filtra as amostras, considerando só aquelas que batem com os dados observados (ABC)
# posterior <- prior[sim_dados == dados_observados] 

## examine o posterior
# hist(posterior) 
# length(posterior)

## verifique se o número de amostras é satisfatório (p. ex., > 1000)

## obtenha alguns sumários do  posterior
# mean(posterior)
# median(posterior)
# sd(posterior)
# mad(posterior)
# quantile(posterior, c(0.025, 0.975))

Questão II) Qual é a probabilidade de que o método A seja melhor que telemarketing?

O Departamento de Marketing acabou de nos dizer que a taxa de inscrição seria de 20% se os assinantes de salmão fossem "capturados" por uma campanha de telemarketing (para nós, não está muito claro de onde o marketing obteve esse número tão preciso). Então, dado o modelo e os dados que obtivemos na última pergunta, qual é a probabilidade de que o método A tenha uma taxa de adesão maior do que o telemarketing?

Dica 1: Se você tem um vetor de amostras representando uma distribuição de probabilidade, que você deve ter da última questão, o cálculo da probabilidade acima de um determinado valor é feito simplesmente contando o número de amostras acima desse valor e dividindo pelo número total de amostras.

Questão III) Se o método A fosse usado em 100 pessoas, qual seria o número de inscrições?

Dica 1: A resposta novamente não é um único número, mas uma distribuição sobre o número provável de inscrições.

Dica 2: Como antes, a distribuição binomial é uma boa candidata para quantas pessoas se inscrevem das 100 possíveis.

Dica 3: Certifique-se de não "jogar fora" a incerteza, por exemplo, usando um resumo da distribuição do posterior calculada na primeira pergunta. Use a amostra original do posterior completa!

Solução (uma entre muitas possíveis):

Questão I

# reprodutibilidade
set.seed(123)

# número de simulações
n_sim <- 10000

# definição do prior de prob e sua amostragem
prior_prob <- runif(n_sim, 0, 1)

# modelo gerativo:
modelo_gerativo <- function(taxa) {
  n_inscritos <- rbinom(1, size = 16, prob = taxa)
  n_inscritos
}

# Simulação dos dados
n_inscritos <- rep(NA, n_sim)
for(i in 1:n_sim) {
  n_inscritos[i] <- modelo_gerativo(prior_prob[i])
}

# Seleção dos valorem de parâmetros que batem com as observações 
post_taxa <- prior_prob[n_inscritos == 6]

# verificando o número se amostras
length(post_taxa)
## [1] 651
# Visualização do posterior.
hist(post_taxa, xlim = c(0, 1),main='',col='gold')
abline(v=mean(post_taxa),lwd=2,col='blue')

# Sumários 
mean(post_taxa)
## [1] 0.3912679
quantile(post_taxa, c(0.025, 0.975))
##      2.5%     97.5% 
## 0.1836081 0.6257492

Questão II

sum(post_taxa > 0.2) / length(post_taxa)
## [1] 0.9539171

Questão III

# Isso pode ser feito com um loop usando for:
assinaturas <- rep(NA, length(post_taxa))
for(i in 1:length(post_taxa)) {
  assinaturas[i] <- rbinom(n = 1, size = 100, prob = post_taxa[i])
}

# Mas, como rbinom é uma função vetorizada, o mesmo procedimento pode ser feito como:
assinaturas <- rbinom(n = length(post_taxa), size = 100, prob = post_taxa)

hist(assinaturas, xlim = c(0, 100),col='gold',main='')
abline(v=mean(assinaturas),lwd=2,col='blue')

quantile(assinaturas, c(0.025, 0.975))
##  2.5% 97.5% 
##    17    63

Em conclusão, uma opinião decente seria que o número de assinaturas ficaria entre 20 e 60.

Vejam também este link: https://www.sumsar.net/blog/2014/10/tiny-data-and-the-socks-of-karl-broman/

Exercícios

Estimativa do parâmetro de uma distribuição exponencial por ABC

Considere os dados

x = c(1.158, 2.220, 0.449, 0.491, 0.758, 0.018, 0.196, 1.461, 3.298, 1.891, 0.583, 0.745, 0.502, 1.065, 0.042, 1.525, 0.651, 0.479, 0.966, 0.795)

Vamos supor que eles obedeçam a uma distribuição exponencial, \(P(x) \propto \exp(-ax)\)

  1. Qual é a verossimilhança dos dados?
  1. Estime \(a\) por máxima verossimilhança.
  1. Estime \(a\) por MCMC assumindo um prior uniforme para \(a\).
  1. Estime \(a\) por MCMC assumindo um prior gaussiano para \(a\), \(N(\mu = 2, \sigma = 2)\).
  1. Compare as medianas das cadeias nos dois casos depois de remover o burn-in (suponha que o burn-in corresponda às 100 primeiras iterações).
  1. Compare também o intervalo de credibilidade de 95% do parâmetro \(a\) depois de remover o burn-in.
  1. Comente como a diferença entre os priores afeta os resultados.
  1. Obtenha um posterior para \(a\) por ABC. Adote um sumário dos dados e um critério de aceitação.