aula de hoje: inferência bayesiana de parâmetros

  1. um problema de astrofísica de altas energias
  1. quantos peixes tem no lago?
  1. modelos gerativos e ABC (Approximate Bayesian computation)
  1. quando os métodos bayesianos e frequentistas divergem?
  1. o modelo beta-binomial

Exercícios

Apêndices

A1. inferência bayesiana: análise sequencial ou conjunta?

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

A3. regressão linear bayesiana

A4. o problema do farol

A5. a estratégia de marketing da Swedish Fish Incorporated

1. 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="contagensX.dat",header=TRUE)
dim(dados)
## [1] 81  2
head(dados)
##     E  N
## 1 2.0 32
## 2 2.1 37
## 3 2.2 18
## 4 2.3 25
## 5 2.4 30
## 6 2.5 22
attach(dados)


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

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.

Em um processo poissoniano, se \(C_i\) é o valor esperado num dado bin, a variância nesse mesmo bin também será \(C_i\), de modo que podemos escrever a verossimilhança como \[L \propto \prod_i \exp \Big [ -{(N_i - C_i)^2 \over C_i} \Big]\] \[log L \propto -\sum_i {(N_i - C_i)^2 \over C_i}\]

Vamos fazer um modelo bayesiano para este processo, assumindo priores uniformes para \(\alpha\) e \(\beta\). Nesse caso o posterior fica proporcional à verossimilhança.

# tempo de processamento
t0 = Sys.time()

# log da verossimilhança dos dados
# param[1] = alfa    param[2] = beta
npar = 2
param=NULL
logver = function(param){
NC=param[1]*E^(-param[2])
#returna valor proporcional ao log da verossimilhança
return(-sum((N-NC)^2/NC))
}

# 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 para o posterior em log
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:
#reprodutibilidade
set.seed(123)
# simulacao
xinicial = c(100,3)
nsim=1000000
simulacao = mcmc(xinicial, nsim)

# vamos supor um burn-in de 10% das iterações:
burnin=0.1*nsim

# taxa de aceitação
aceitacao = 1-mean(duplicated(simulacao[-(1:burnin),]))
aceitacao
## [1] 0.3044174
#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. 
##   131.6   174.5   186.1   186.8   198.5   254.8
summary(p2)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.112   2.373   2.424   2.423   2.475   2.720
# intervalo de de credibilidade de 95%: a probabilidade do parâmetro estar dentro deste intervalo é 95%
quantile(p1, probs = c(0.025, 0.975))
##     2.5%    97.5% 
## 153.4755 223.5588
quantile(p2, probs = c(0.025, 0.975))
##     2.5%    97.5% 
## 2.272523 2.567655
# desvio padrão
c(sd(p1),sd(p2))
## [1] 17.76428231  0.07520206
# moda
estimate_mode <- function(x) {
  d <- density(x,bw='sj')
  d$x[which.max(d$y)]
}
m1 = estimate_mode(p1)
m2 = estimate_mode(p2)

# moda dos parâmetros
c(m1,m2)
## [1] 184.54999   2.42416
# 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')

# tempo de processamento
Sys.time()-t0
## Time difference of 43.00833 secs

2. 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, 7
  • \(K\) o número de peixes pescados, 20
  • \(X\) o número dentre eles que estão marcados, 3

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

# tempo de processamento
t0 = Sys.time()

# reprodutibilidade
set.seed(123)

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

# solução de máxima verossimilhança
MV = M*K/X
MV
## [1] 46.66667
# queremos T, o número total de peixes no lago
# 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',lwd=2,lty=3)
abline(v=T975,col='blue',lwd=2,lty=3)

# moda:
T[y == max(y)]
## [1] 46
# tempo de processamento
Sys.time()-t0
## Time difference of 0.01992559 secs

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

3. modelos gerativos e ABC (Approximate Bayesian Computation)

modelos gerativos: inferência baseada em simulações dos dados:

  • simulamos parâmetros a partir de um prior, \(P(w)\)
  • usamos esses parâmetros para simular “dados”, e
  • usamos os dados simulados para estimar o posterior dos parâmetros

Vamos agora resolver o problema dos peixes no lago com modelos gerativos e usando ABC. Note que não vamos precisar calcular a verossimilhança!

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

# tempo de processamento
t0 = Sys.time()

# reprodutibilidade
set.seed(123)

# dados
# número de peixes marcados
M=7
# número de peixes pescados
K=20
# número de peixes pescados marcados
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%
q = quantile(post_npeixes, c(0.025, 0.975))
q
##  2.5% 97.5% 
##    31   298
abline(v=q[[1]],col='blue',lwd=2, lty=3)
abline(v=q[[2]],col='blue',lwd=2, lty=3)

# compare com o intervalo obtido acima

#tempo de processamento
Sys.time()-t0
## Time difference of 0.9396336 secs

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

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

5. 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, \(p\): a probabilidade de sucesso em cada tentativa (\(0 \le p \le 1\))

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

Dado um conjunto de dados, \(D=\{n,N\}\), vamos fazer um modelo bayesiano para estimar \(p\)

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

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

O posterior então fica: \[ P(p|D) \propto p^{n+a-1}(1-p)^{N-n+b-1} \] ou \[ P(p|D) = Beta(p|n+a,N-n+b) \]

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

observação: a função beta e a distribuição beta são coisas diferentes!

  • função beta: \[\beta(a,b) = \int_0^1 p^{a-1}(1-p)^{b-1}dp = \frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}\]
  • distribuição beta: \[P(p|a,b) = \frac{1}{\beta(a,b)}p^{a-1}(1-p)^{b-1}\]

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 são 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 \(p\) 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:

# tempo de processamento
t0 = Sys.time()

# reprodutibilidade
set.seed(123)

# número de galáxias na amostra
N = 1000

# número de ET na amostra
n = 735

# probabilidade de ET na amostra
pET = n/N

# sequência com o número de galáxias 
x = seq(1, N, 1)

# verossimilhança
ver = dbinom(x, size = N, prob = pET)

# visualização
plot(x, ver, type = "l", lwd = 3, main = "verossimilhança", ylab = "densidade de probabilidades", xlab = "número de galáxias early type", col='blue')
abline(v = n, col = "red",lwd=2)

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

a=1
b=1
x=seq(0,1,0.001)
prior = dbeta(x,a,b)

plot(x,prior,type='l',col='blue', main='prior', ylab = "densidade de probabilidades",xlab=expression(p),lwd=3)

Nesse caso, o posterior será: \[ P(p|D,a,b) =Beta(p|n+a,N-n+b) = Beta(p|a=736,b=266)\]

post = dbeta(x,n+a,N-n+b)

plot(x,post,type='l',col='blue', main='posterior', ylab = "densidade de probabilidades",xlab=expression(p),lwd=3)

# moda  
moda = x[which(post == max(post))]
moda
## [1] 0.735
# intervalo de de credibilidade de 95%:
# para isso precisamos determinar os valores de p correspondentes aos quantis 0.025 e 0.975

# defino uma função de distribuição de probabilidades para o posterior
f = function(p) dbeta(p,n+a,N-n+b)

# defino a função cumulativa
cum_prob = function(x) integrate(f,0,x)$value

# defino uma função para o quantil inferior
h = function(x) cum_prob(x) - 0.025

# resolvo a equação h=0 para determinar o quantil inferior
q1 = uniroot(h, interval=c(0, 1))$root

# repito o procedimento para o quantil superior
h = function(x) cum_prob(x) - 0.975
q2 = uniroot(h, interval=c(0, 1))$root

# resultado
round(c(q1,q2),3)
## [1] 0.707 0.761

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 (\(a\) e \(b\), os “shape parameters”) que podem ser obtidos da média e variância da distribuição:

estBetaParams <- function(media, variancia) {
  a <- ((1 - media) / variancia - 1 / media) * media^2
  b <- a * (1 / media - 1)
  return(c(a,b))
}

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

prior.params <- estBetaParams(media = 0.6, variancia = 0.2^2)
prior.params
## [1] 3 2

Vamos representar graficamente este prior:

x <- seq(from = 0,to = 1,  by = 0.001)
# note que seq deve estar no intervalo da função Beta, [0,1]

# prior
prior = dbeta(x, shape1 = prior.params[[1]], shape2 = prior.params[[2]])
plot(x, prior, lwd = 3, col = 'blue', ylab = "densidade", xlab = expression(p), type = "l", main = "prior")

e o posterior será:

N <- 1000
n <- 735

# hiper-parâmetros do posterior
a_post <- prior.params[[1]] + n
b_post <- prior.params[[2]] + N - n

# posterior
post = dbeta(x, shape1 = a_post, shape2 = b_post)

plot(x, post, lwd = 3, type = "l", ylab = "densidade", xlab = expression(p), col = "red", main = "posterior e prior")  

# vamos plotar o prior junto com o posterior
lines(x, prior, 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
moda = x[which(post == max(post))]
moda
## [1] 0.735
# intervalo de credibilidade de 95%:
# aplicamos o mesmo procedimento descrito acima
f = function(x) dbeta(x, shape1 = a_post, shape2 = b_post)
cum_prob = function(x) integrate(f,0,x)$value
h = function(x) cum_prob(x) - 0.025
q1 = uniroot(h, interval=c(0, 1))$root
h = function(x) cum_prob(x) - 0.975
q2 = uniroot(h, interval=c(0, 1))$root
round(c(q1,q2),3)
## [1] 0.707 0.761
# tempo de processamento
Sys.time()-t0
## Time difference of 0.4622991 secs

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 \(p\) nos dois casos. Se o número de dados fosse bem menor as diferenças seriam maiores!

Exercícios

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

Inicialize a semente de números aleatórios com seu número USP

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)\), \(a > 0\)

  1. Qual é a verossimilhança dos dados?
  1. Estime \(a\) por máxima verossimilhança.
  1. Estime o posterior de \(a\) por MCMC, assumindo um prior uniforme para \(a\). Lembre-se que a função proposta deve fornecer apenas \(a>0\).
  1. Estime o posterior de \(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.

Apêndice

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

  • temos um conjunto de dados \(D=\{x_i,\sigma\}\) 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, já que o posterior, a distribuição de probabilidades do parâmetro \(\mu\), é \[P(\mu|D) \propto P(D|\mu) P(\mu)\]
  • verossimilhança dos dados: \[ P(D|\mu) = \prod_{i=1}^n P(x_i,\sigma|\mu) \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: \(P(\mu|D) = N( \mu_n, \tau_n)\), onde \[ \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 migra 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

A2. 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 a partir desses dados.

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.

# tempo de processamento
t0 = Sys.time()

# 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])
dados
##  [1] 1.4395244 1.7698225 3.5587083 2.0705084 2.1292877 3.7150650 2.4609162
##  [8] 0.7349388 1.3131471 1.5543380
# 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
# vamos adotar:
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)

t0 = Sys.time()
posterior <- MCMC_ABC(parini,niter)
Sys.time() - t0
## Time difference of 12.02932 secs
medias = posterior[,1]
desvp = posterior[,2]

# moda
m1 = estimate_mode(posterior[,1])
m2 = estimate_mode(posterior[,2])
c(m1,m2)
## [1] 2.0896849 0.9619098
# como estamos fazendo muitas simulações o burn-in não é importante

# valores médios dos parâmetros simulados
c(mean(medias),mean(desvp))
## [1] 2.062755 1.112628
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
# tempo de processamento
Sys.time()-t0
## Time difference of 12.92981 secs

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

A3. regressão linear bayesiana

Vamos fazer uma regressão linear bayesiana.

Vamos gerar dados a partir de um modelo \(y = 1 + 2x\) com erros gaussianos de desvio padrão 0.5

# tempo de processamento
t0 = Sys.time()

# modelo y=a+bx+e
a = 1
b = 2
sigma = 0.5

# simulação dos dados:
set.seed(2903)
n = 10
x = runif(n,min=0,max=5)
y = a+b*x+rnorm(n,mean=0,sd=sigma)

solução de máxima verossimilhança:

# máxima verossimilhança
mv = lm(y~x)
mv
## 
## Call:
## lm(formula = y ~ x)
## 
## Coefficients:
## (Intercept)            x  
##       1.424        1.839
# parâmetros
aMV = mv$coeff[1]
bMV = mv$coeff[2]

# visualização
plot(x,y,pch=20, col='salmon',main='M: y = 1 + 2x + eps',cex=2)
abline(a=1,b=2,lwd=2, col='green')
abline(a=mv$coeff[1],b=mv$coeff[2],lwd=2,col='red')
legend('topleft',legend=c('modelo','ajuste bayesiano'),col=c('green','red'),  text.col = c('green','red'),bty='n')
legend('bottomright',legend=c('MV:     a=0.7268  b=2.0895  sigma=0.4504','bayes: a=0.7160  b=2.0883  sigma=0.3820'),col=c('blue','red'),  text.col = c('blue','red'),bty='n')

Para a solução bayesiana, vamos adotar, priores uniformes para \(a\) e \(b\) e o prior de Jeffreys (uniforme em log) para \(\sigma\): \[P(a,b,\sigma) \propto {1 \over \sigma^2}\]

Nesse caso o posterior fica \[P(a,b,\sigma|D) \propto \sigma^{-(n+2)}\exp \Big[- \sum_i{(y_i - a-b x_i)^2 \over 2 \sigma^2} \Big]\]

Vamos resolver com MCMC

# log do posterior
# param[1] = a, param[2] = b, param[3] = sigma
npar = 3
param = rep(0,npar)
logposterior = function(param){
del = y - param[1]-param[2]*x
return(-(n+2)*log(param[3])-1/(2*param[3]^2)*sum(del^2))
}

# função proposta gaussiana para cada parâmetro
 pp = param
funcao_proposta = function(param){
pp[1] = rnorm(1,mean = param[1], sd= 0.1)
pp[2] = rnorm(1,mean = param[2], sd= 0.1)
# para sigma: amostragem gaussiana em torno de log(sigma)
pp[3] = exp(rnorm(1,mean = log(param[3]), sd= 0.1))
return(pp)
}

# mcmc 
# inicialização dos parâmetros
xini = c(1,1,1)
# tamanho da cadeia
niter = 1000000
# mcmc
set.seed(234234)
cadeiamc = mcmc(xini, niter)

# burn-in
burnin = 0.1*niter
aceitacao = 1-mean(duplicated(cadeiamc[-(1:burnin),]))
aceitacao
## [1] 0.3288907
# visualização sem o burn-in
# notem como se pode remover o burn-in da cadeia
par(mfrow = c(2,3))
hist(cadeiamc[-(1:burnin),1],nclass=30, main="distribuição de a", xlab='', ylab="frequência",col='gold',cex.lab=1.2)
abline(v = mean(cadeiamc[-(1:burnin),1]),lwd=2)
hist(cadeiamc[-(1:burnin),2],nclass=30, main="distribuição de b", xlab='', ylab="frequência",col='gold',cex.lab=1.2)
abline(v = mean(cadeiamc[-(1:burnin),2]),lwd=2)
hist(cadeiamc[-(1:burnin),3],nclass=30, main="distribuição de sigma", xlab='', ylab="frequência",col='gold',cex.lab=1.2)
abline(v = mean(cadeiamc[-(1:burnin),3]),lwd=2)

plot(cadeiamc[-(1:burnin),1], type = "l", xlab="Cadeia de  valores de a" , main = "", col='gray27',cex.lab=1.2,ylab='a')
abline(h = mean(cadeiamc[-(1:burnin),1]), col="red",lwd=2)

plot(cadeiamc[-(1:burnin),2], type = "l", xlab="Cadeia de  valores de b" , main = "", col='gray27',cex.lab=1.2,ylab='b')
abline(h = mean(cadeiamc[-(1:burnin),2]), col="red",lwd=2)

plot(cadeiamc[-(1:burnin),1], type = "l", xlab="Cadeia de  valores de sigma" , main = "", col='gray27',cex.lab=1.2,ylab='sigma')
abline(h = mean(cadeiamc[-(1:burnin),3]), col="red",lwd=2)

# estatísticas
# valores médios:
mean(cadeiamc[-(1:burnin),1])
## [1] 1.423363
mean(cadeiamc[-(1:burnin),2])
## [1] 1.839443
mean(cadeiamc[-(1:burnin),3])
## [1] 0.3435634
# valores medianos:
median(cadeiamc[-(1:burnin),1])
## [1] 1.423515
median(cadeiamc[-(1:burnin),2])
## [1] 1.839393
median(cadeiamc[-(1:burnin),3])
## [1] 0.3285368
# moda
m1 = estimate_mode(cadeiamc[-(1:burnin),1])
m2 = estimate_mode(cadeiamc[-(1:burnin),2])
m3 = estimate_mode(cadeiamc[-(1:burnin),3])
c(m1,m2,m3)
## [1] 1.4319361 1.8457029 0.3059936
# intervalo de credibilidade dos parâmetros
# 90%
# tem 90% de probabilidade de que o valor verdadeiro esteja dentro desse intervalo
quantile(cadeiamc[-(1:burnin),1], probs=c(0.05,0.95))
##       5%      95% 
## 1.069949 1.777287
quantile(cadeiamc[-(1:burnin),2], probs=c(0.05,0.95))
##       5%      95% 
## 1.724569 1.954374
quantile(cadeiamc[-(1:burnin),3], probs=c(0.05,0.95))
##        5%       95% 
## 0.2348184 0.5028114
# distribuições conjuntas
par(mfrow = c(1,3))
smoothScatter(cadeiamc[-(1:burnin),1],cadeiamc[-(1:burnin),2],xlab='a',ylab='b',cex.lab=1.5)
smoothScatter(cadeiamc[-(1:burnin),1],cadeiamc[-(1:burnin),3],xlab='a',ylab='sigma',cex.lab=1.5)
smoothScatter(cadeiamc[-(1:burnin),2],cadeiamc[-(1:burnin),3],xlab='b',ylab='sigma',cex.lab=1.5)

# resultado
par(mfrow = c(1,1))
plot(x,y,pch=20, col='salmon',main='M: y = 1 +2x + eps',cex=2)
abline(a=1,b=2,lwd=2, col='green')
abline(a=mv$coeff[1],b=mv$coeff[2],lwd=2)
abline(a=mean(cadeiamc[-(1:burnin),1]),b=mean(cadeiamc[-(1:burnin),2]),lwd=2, col='blue')

ntot=100000
xx = rep(0,ntot)
yy = rep(0,ntot)

for(i in 1:ntot){
k = sample(seq(burnin+1,niter,1),1)
xx[i]=runif(1,min=0,max=5)
yy[i]=cadeiamc[k,1]+cadeiamc[k,2]*xx[i]
}

d = yy - mean(cadeiamc[-(1:burnin),1])-mean(cadeiamc[-(1:burnin),2])*xx
q = quantile(d,probs=c(0.005,0.05,0.25,0.75,0.95,0.995))
hist(d,breaks=100)
abline(v=q,lty=3,col='blue')

smoothScatter(xx,yy,main='M: y = 1 + 2x + eps',xlab='x',ylab='y')
points(x,y, col='salmon',main='M: y = 1 +2x + eps',cex=2,pch=20)
abline(a=1,b=2,lwd=2, col='green')
abline(a=mv$coeff[1],b=mv$coeff[2],lwd=2)
abline(a=mean(cadeiamc[-(1:burnin),1]),b=mean(cadeiamc[-(1:burnin),2]),lwd=2, col='red')
legend('topleft',legend=c('modelo','ajuste bayesiano'),col=c('green','red'),  text.col = c('green','red'),bty='n')
for(i in 1:length(q)){
segments(0,mean(cadeiamc[-(1:burnin),1])+q[i],5,mean(cadeiamc[-(1:burnin),1])+5*mean(cadeiamc[-(1:burnin),2])+q[i])
}
segments(3,2,4,2)
legend(2,2,legend = '0.5%,5%,25%,75%,95%,99.5%',bty='n')

# tempo de processamento
Sys.time() - t0
## Time difference of 9.604361 mins

A4. 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 = \{x_1, x_2, ..., x_n\}\), 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; para um dado \(x_i\), \[P(x_i|\alpha,\beta) = {\beta \over \pi [\beta^2 + (x_i-\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:

# tempo de processamento
t0 = Sys.time()

# 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 p 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: como as medidas são independentes, o log da verossimilhança é a soma dos logs das verossimilhanças de cada dado individualmente \[{\cal L}(\alpha,\beta) = \prod_{k=1}^N P(x_k|\alpha,\beta) = \prod_{k=1}^N {\beta \over \pi [\beta^2 + (x_i-\alpha)^2]}\] e \[\ln {\cal L}(\alpha,\beta) = \sum_{k=1}^N\ln {\beta \over \pi [\beta^2 + (x_i-\alpha)^2]}\]

# 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: vamos supor que \(\alpha\) e \(\beta\) são independentes e com priores uniformes \[P(\alpha,\beta) = P(\alpha)P(\beta) = {\rm cte}\] em log: \[\ln P(\alpha,\beta) = \ln P(\alpha) + \ln P(\beta)\]

# 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)}

Vamos amostrar este posterior com MCMC.

Para isso vamos adotar uma 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 (ver acima):

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

# vamos supor um burn-in de 10% das 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=0.1*nsim

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

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.4429398
media.b=mean(cb)
media.b
## [1] 1.514905
# desvios padrão
sd(ca)
## [1] 0.2160158
sd(cb)
## [1] 0.2138667
# intervalo de credibilidade dos parâmetros:
quantile(ca, c(0.025, 0.975))
##       2.5%      97.5% 
## 0.02109928 0.87124807
quantile(cb, c(0.025, 0.975))
##     2.5%    97.5% 
## 1.136961 1.974206

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)
## Warning in if (!add) {: the condition has length > 1 and only the first element
## will be used
## Warning in if (!add) {: the condition has length > 1 and only the first element
## will be used

# vermelho: parâmetros médios da simulação       verde: parâmetros verdadeiros, usados para gerar os dados

# tempo de processamento
Sys.time()-t0
## Time difference of 2.00867 mins

A5. 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 peixes 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:

# tempo de processamento
t0 = Sys.time()
## 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
# tempo de processamento
Sys.time()-t0
## Time difference of 0.1669312 secs

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/