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

  1. um problema de astrofísica de altas energias
  1. regressão linear bayesiana
  1. quantos peixes tem no lago?
  1. quando os métodos bayesianos e frequentistas divergem…
  1. o modelo beta-binomial

Exercícios

Apêndice: o problema do farol

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
colnames(dados)
## [1] "E" "N"
# vamos acessar as variáveis contidas em "dados" pelo nome
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\):

\[20 < \alpha < 300, ~~~~1.5 < \beta < 3\]

Nesse caso o posterior fica proporcional à verossimilhança.

Vamos estimar o posterior por MCMC:

# tempo de processamento deste módulo
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
# a cruz vermelha corresponde ao valor médio dos parâmetros
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 39.03851 secs

Note que todo modelo envolve escolhas! No caso, poderíamos, por exemplo, analisar \(log N \times log E\), que é um problema linear,

ou poderíamos supor os erros não conhecidos,

ou …

George Box: “all models are wrong, but some are useful”

2. regressão linear bayesiana

Aparentemente toda galáxia massiva abriga em seu núcleo um buraco negro supermassivo (SMBH) cuja massa se relaciona com certas propriedades globais das galáxias.

McConnell & Ma (2013, ApJ, 764, 184) fizeram uma compilação de dados de 72 SMBHs e de suas galáxias hospedeiras para estudar relações de escala entre a massa do SMBH e propriedades das galáxias.

Vamos aqui analisar a relação entre a massa do SMBH e a dispersão central de velocidades \(\sigma\) da galáxia.

\(\sigma\) é uma medida da profundidade do poço de potencial da galáxia e, no caso das galáxias elípticas, pode ser usada para estimar a massa total (relação de Faber-Jackson).

leitura dos dados:

# tempo de processamento deste módulo
t0 = Sys.time()

dados = read.table(file='McConnell_Ma_2013_ascii.txt',head = TRUE)
dim(dados)
## [1] 73 15
colnames(dados)
##  [1] "Galaxy"     "DMpc"       "MBH"        "MBHd"       "MBHu"      
##  [6] "method"     "sigma"      "sigmad"     "sigmau"     "logLV"     
## [11] "elogLV"     "Mbulge"     "rinf"       "morphology" "profile"
# descrição das colunas:
# "Galaxy"    nome da galáxia
# "DMpc"      distância (Mpc)
# "MBH"       massa do BN 
# "MBHd"      massa inferior do BN (68%) 
# "MBHu"      massa superior do BN (68%)
# "method"    observável usado para estimar a massa do BN 
# "sigma"     dispersão de velocidades das estrelas (km/s) 
# "sigmad"    valor inferior para sigma (68%)
# "sigmau"    valor superior para sigma (68%)
# "logLV"     log da luminosidade da galáxia na banda V (em unidades de Lsun)
# "elogLV"    erro em logLV 
# "Mbulge"    massa do bojo da galáxia (em Msun)
# "rinf"      "raio de influência" do BN (arcsec)  GM_BN/sigma^2
# "morphology" tipo morfológico
# "profile"   tipo de perfil central: "core" (C), lei de potência (pl)

# note que os erros em MBH e sigma são assimétricos!

attach(dados)
n = nrow(dados)

# dados em log
# x: log sigma
# y: log MBH
x = log10(sigma)
xd = log10(sigma)-log10(sigmad) # barra de erro em x inferior
xu = log10(sigmau)-log10(sigma) # barra de erro em x superior
y = log10(MBH)
yd = log10(MBH)-log10(MBHd) # barra de erro em y inferior
yu = log10(MBHu)-log10(MBH) # barra de erro em x superior

# função para visualização com barras de erro em x e y
# DeepSeek
#' Scatter Plot with Asymmetrical Error Bars
#'
#' @param x Vector of x values
#' @param y Vector of y values
#' @param x_lower Vector of lower x errors
#' @param x_upper Vector of upper x errors
#' @param y_lower Vector of lower y errors
#' @param y_upper Vector of upper y errors
#' @param xlab X-axis label (default: "X")
#' @param ylab Y-axis label (default: "Y")
#' @param main Plot title (default: "")
#' @param pch Point character (default: 19)
#' @param col Point color (default: "black")
#' @param err_col Error bar color (default: "gray50")
#' @param err_lwd Error bar line width (default: 1)
#' @param err_length Error bar cap length (in inches, default: 0.03)
#' @param ... Additional graphical parameters passed to plot()
#' @return A scatter plot with asymmetrical error bars
#' @examples
#' # Example usage:
#' x <- 1:5
#' y <- c(2, 4, 3, 5, 6)
#' x_lower <- c(0.1, 0.2, 0.1, 0.3, 0.2)
#' x_upper <- c(0.2, 0.3, 0.2, 0.4, 0.3)
#' y_lower <- c(0.3, 0.4, 0.2, 0.5, 0.4)
#' y_upper <- c(0.4, 0.5, 0.3, 0.6, 0.5)
#' plot_with_asym_errors(x, y, x_lower, x_upper, y_lower, y_upper)
plot_with_asym_errors <- function(x, y, 
                                 x_lower = NULL, x_upper = NULL, 
                                 y_lower = NULL, y_upper = NULL,
                                 xlab = "X", ylab = "Y", main = "",
                                 pch = 19, col = "black",
                                 err_col = "gray50", err_lwd = 1,
                                 err_length = 0.03, ...) {
  
  # Input validation
  if (length(x) != length(y)) {
    stop("x and y must have the same length")
  }
  
  # Set default errors if not provided
  if (is.null(x_lower)) x_lower <- rep(0, length(x))
  if (is.null(x_upper)) x_upper <- rep(0, length(x))
  if (is.null(y_lower)) y_lower <- rep(0, length(x))
  if (is.null(y_upper)) y_upper <- rep(0, length(x))
  
  # Check error lengths
  if (length(x_lower) != length(x) || length(x_upper) != length(x) ||
      length(y_lower) != length(x) || length(y_upper) != length(x)) {
    stop("All error vectors must have same length as x and y")
  }
  
  # Calculate plot limits
  x_range <- range(c(x - x_lower, x + x_upper), na.rm = TRUE)
  y_range <- range(c(y - y_lower, y + y_upper), na.rm = TRUE)
  
  # Create base plot
  plot(x, y, xlim = x_range, ylim = y_range, 
       xlab = xlab, ylab = ylab, main = main,
       pch = pch, col = col, ...)
  
  # Add error bars
  arrows(x0 = x - x_lower, y0 = y,
         x1 = x + x_upper, y1 = y,
         angle = 90, code = 3, length = err_length,
         col = err_col, lwd = err_lwd)
  
  arrows(x0 = x, y0 = y - y_lower,
         x1 = x, y1 = y + y_upper,
         angle = 90, code = 3, length = err_length,
         col = err_col, lwd = err_lwd)
  
  # Return invisibly
  invisible(list(x = x, y = y, 
                x_err = cbind(x_lower, x_upper),
                y_err = cbind(y_lower, y_upper)))
}

par(mfrow = c(1,1))
library(latex2exp)
plot_with_asym_errors(x, y, xd, xu, yd, yu,
                     xlab = expression(paste('log ',sigma,' (km/s)')), ylab = TeX('$log (M_{BH} /M_{sun})$'),
                     main = " ",
                     pch = 21, col = "blue", bg = "lightblue",
                     err_col = "lightblue", err_lwd = 1.2)

Vamos considerar diversos modelos e comparar os resultados.

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

# a função lm() determina a solução de máxima verossimilhança para o ajuste linear
# vamos ignorar os erros
mv = lm(y~x)

# solução
summary(mv)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.51463 -0.29059  0.06777  0.24989  0.95344 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -4.4266     0.7280   -6.08 5.43e-08 ***
## x             5.5359     0.3167   17.48  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4377 on 71 degrees of freedom
## Multiple R-squared:  0.8114, Adjusted R-squared:  0.8088 
## F-statistic: 305.5 on 1 and 71 DF,  p-value: < 2.2e-16
# parâmetros
aMV = mv$coeff[1]
bMV = mv$coeff[2]
sMV = sd(mv$residuals)

# visualização
plot(x,y,pch=20, col='lightblue',main='máxima verossimilhança',cex=2)
abline(a=mv$coeff[1],b=mv$coeff[2],lwd=2,col='green')

Modelo 2: uma solução bayesiana

modelo com parâmetros \(a,b,\sigma\) (\(\sigma\): incerteza nas estimativas de \(y\)) e com priores uniformes para \(a\) e \(b\) e o prior de Jeffreys para \(\sigma\): \[P(a,b,\sigma) \propto {1 \over \sigma^2}\]

Nesse caso o posterior fica (\(n\): número de pontos) \[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 amostrar o posterior 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.1857998
# 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),3], 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:
aMB1 = mean(cadeiamc[-(1:burnin),1])
bMB1 = mean(cadeiamc[-(1:burnin),2])
sMB1 = mean(cadeiamc[-(1:burnin),3])

# 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)
points(aMB1,bMB1,pch=3,col="red",lwd=3)   
smoothScatter(cadeiamc[-(1:burnin),1],cadeiamc[-(1:burnin),3],xlab='a',ylab='sigma',cex.lab=1.5)
points(aMB1,sMB1,pch=3,col="red",lwd=3) 
smoothScatter(cadeiamc[-(1:burnin),2],cadeiamc[-(1:burnin),3],xlab='b',ylab='sigma',cex.lab=1.5)
points(bMB1,sMB1,pch=3,col="red",lwd=3) 

# intervalo de credibilidade dos parâmetros
# 68%
# tem 68% de probabilidade de que o valor verdadeiro esteja dentro desse intervalo
qaMB1 = quantile(cadeiamc[-(1:burnin),1], probs=c(0.16,0.84))
qbMB1 = quantile(cadeiamc[-(1:burnin),2], probs=c(0.16,0.84))
qsMB1 = quantile(cadeiamc[-(1:burnin),3], probs=c(0.16,0.84))

# para comparar com os erros da MV
sa = as.numeric((qaMB1[2]-qaMB1[1])/2)
sb = as.numeric((qbMB1[2]-qbMB1[1])/2)
ss = as.numeric((qsMB1[2]-qsMB1[1])/2)

# parâmetro     valor    erro
print(c('a: ',round(aMB1,3),round(sa,3)))
## [1] "a: "    "-4.436" "0.71"
print(c('b: ',round(bMB1,3),round(sb,3)))
## [1] "b: "   "5.54"  "0.309"
print(c('sigma: ',round(sMB1,3),round(ss,3)))
## [1] "sigma: " "0.436"   "0.036"
# visualização
par(mfrow = c(1,1))
plot(x,y,pch=20, col='lightblue',main='MB1',cex=2)
abline(aMV,b=bMV,lwd=2,col='green')
abline(aMB1,bMB1,lwd=2, col='blue')

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

Modelo 3: modelo gerativo para erros em x e y

Vamos usar um modelo gerativo para

  • simularmos dados dada uma verossimilhança
  • para cada conjunto de dados simulados estimamos os parâmetros da regressão linear \((a,b,\sigma)\)
  • usamos a distribuição de \((a,b,\sigma)\) para estimar o posterior destes parâmetros

Vamos supor que um “novo” dado possa ser gerado de um ponto \((x,y)\) como \[P(x) \sim P(x_{obs},\sigma_{x1},\sigma_{x2})~~~~~~~~~~~P(y) \sim P(y_{obs},\sigma_{y1},\sigma_{y2})\] que são as verossimilhanças de cada \(x\) e \(y\) dados os valores observados \((x_{obs},y_{obs})\) e seus erros (assimétricos!)

Vamos escrever a verossimilhança de \(x\), com erros assimétricos, como \(L(x \mid x_{obs})\): \[ L(x \mid x_{obs}) = \begin{cases} \displaystyle \frac{1}{\sqrt{2\pi}\,\sigma_1} \exp\!\left(-\frac{(x_{obs}-x)^2}{2\sigma_1^2}\right), & \text{se } x_{obs}< x, \\[1ex] \displaystyle \frac{1}{\sqrt{2\pi}\,\sigma_2} \exp\!\left(-\frac{(x_{obs}-x)^2}{2\sigma_2^2}\right), & \text{se } x_{obs}\ge x, \end{cases} \] onde \(\sigma_1\) é a incerteza (o desvio padrão) quando \(x\) cai abaixo de \(x_{obs}\) e \(\sigma_2\) quando \(x\) é maior ou igual a \(x_{obs}\)

Idem para \(L(y \mid y_{obs})\)

# tempo de processamento
t0 = Sys.time()

# reprodutibilidade
set.seed(123)

# amostragem de um ponto para um dado com barras de erro assimétricas
#D_obs: valor observado, sigma1: erro para D < D_obs, sigma2: erro para D >= D_obs
amostragem_erros_assim <- function(D_obs, sigma1, sigma2) {
# decidimos de que lado vamos amostrar  
  lado = runif(1)
  amostra <- ifelse(lado < 0.5,
                       rnorm(1, mean = D_obs, sd = sigma1),
                       rnorm(1, mean = D_obs, sd = sigma2))
  return(amostra)
}


# para cada um dos n pontos faremos nsim simulações
# vetores para armazernar resultados
nsim = 10000
xs = rep(0,n)
ys = rep(0,n)
as = rep(0,nsim)
bs = rep(0,nsim)
sigmas = rep(0,nsim)

# para cada simulação
for(j in 1:nsim){
# amostramos n valores (xx,ys)  
  for(i in 1:n){
  xs[i] = amostragem_erros_assim(x[i], xd[i], xu[i])
  ys[i] = amostragem_erros_assim(y[i], yd[i], yu[i])
  }
# ajustamos uma reta a esses valores e obtemos os parâmetros  
  ols = lm(ys ~ xs)
  as[j] <- ols$coef[1]
  bs[j] <- ols$coef[2]
  sigmas[j] <- sd(ols$residuals)
}

# distribuição dos parâmetros
par(mfrow = c(1,3))
hist(as,xlab='a',main='')
hist(bs,xlab='b',main='')
hist(sigmas,xlab='sigma',main='')

smoothScatter(as,bs,xlab='a',ylab='b',cex.lab=1.5)
points(mean(as),mean(bs),pch=3,col="red",lwd=3)   
smoothScatter(as,sigmas,xlab='a',ylab='sigma',cex.lab=1.5)
points(mean(as),mean(sigmas),pch=3,col="red",lwd=3)   
smoothScatter(bs,sigmas,xlab='b',ylab='sigma',cex.lab=1.5)
points(mean(bs),mean(sigmas),pch=3,col="red",lwd=3)   

# estatísticas
# valores médios:
aMG = mean(as)
bMG = mean(bs)
sMG = mean(sigmas)

# intervalo de credibilidade dos parâmetros
# 68%
# tem 68% de probabilidade de que o valor verdadeiro esteja dentro desse intervalo
qaMG = quantile(as, probs=c(0.16,0.84))
qbMG = quantile(bs, probs=c(0.16,0.84))
qsMG = quantile(sigmas, probs=c(0.16,0.84))

# para comparar com os erros da MV
sa = as.numeric((qaMG[2]-qaMG[1])/2)
sb = as.numeric((qbMG[2]-qbMG[1])/2)
ss = as.numeric((qsMG[2]-qsMG[1])/2)

# parâmetro     valor    erro
print(c('a: ',round(aMG,3),round(sa,3)))
## [1] "a: "    "-3.968" "0.406"
print(c('b: ',round(bMG,3),round(sb,3)))
## [1] "b: "   "5.336" "0.174"
print(c('sigma: ',round(sMG,3),round(ss,3)))
## [1] "sigma: " "0.485"   "0.027"
# visualização
par(mfrow = c(1,1))
plot(x,y,pch=20, col='lightblue',main='MG',cex=2)
abline(aMV,b=bMV,lwd=2,col='green')
abline(aMB1,bMB1,lwd=2, col='blue')
abline(aMG,bMG,lwd=2, col='salmon')
legend('topleft',legend=c('MV','MB1','MG'),text.col=c('green','blue','salmon'),  bty='n')

# resumo dos resultados com os vários modelos
c(aMV,aMB1,aMG)
## (Intercept)                         
##   -4.426560   -4.436224   -3.967811
c(bMV,bMB1,bMG)
##        x                   
## 5.535944 5.540161 5.335668
c(sMV,sMB1,sMG)
## [1] 0.4346413 0.4356914 0.4852848
# tempo de processamento
Sys.time() - t0
## Time difference of 11.74146 secs

3. 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.018785 secs

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

Vamos agora resolver o problema dos peixes no lago com um modelo gerativo 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)

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
# 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.9951231 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.2852306 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.4915345 secs

Pode-se notar que, nesse exemplo, a forma do prior afeta muito pouco o posterior. Compare a moda e o intervalo de credibilidade 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 (N)

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,\sigma)\). Faça \(set.seed=N\) e em seguida \(\mu = 2 + 0.5*(runif(1)-0.5), ~~ \sigma = 2+runif(1)\).
  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: 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")  

# 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.0896 mins