aula 6

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

1. inferência sequencial ou conjunta?

2. modelos conjugados

3. um problema de astrofísica de altas energias

4. regressão linear bayesiana

5. o problema do farol

6. quantos peixes tem no lago?

7. Exercício

---

---

1. inferência sequencial ou conjunta (bayesian update)?

Vamos ilustrar a inferência sequencial usando a distribuição 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) \]

Para visualizar a inferência sequencial, vamos usar um código criado por Rasmus Bååth e disponível em
https://www.sumsar.net/blog/2018/12/visualizing-the-beta-binomial/

A função abaixo recebe um número de sucessos e falhas codificados como um vetor VERDADEIRO/FALSO ou 0/1. O resultado é uma visualização de como um modelo Beta-Binomial aprende gradualmente a proporção subjacente de sucessos usando esses dados.

A função também retorna uma amostra do posterior, que pode ser usada para, por exemplo, estimar intervalos de credibilidade do parâmetro \(p\).

Vamos usar um prior uniforme entre 0 e 1, que pode ser representado por uma distribuição Beta(1,1); isso pode ser modificado usando-se o argumento prior_prop.

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}\]

# reprodutibilidade
set.seed(123)

# Example usage:
# data <- c(TRUE, FALSE, TRUE, TRUE, FALSE, TRUE, TRUE)
# prop_model(data)
prop_model <- function(data = c(), prior_prop = c(1, 1), n_draws = 10000) {
  suppressPackageStartupMessages(library(tidyverse))
  data <- as.logical(data)
  # data_indices decides what densities to plot between the prior and the posterior
  # For 20 datapoints and less we're plotting all of them.
  data_indices <- round(seq(0, length(data), length.out = min(length(data) + 1, 20)))

  # dens_curves will be a data frame with the x & y coordinates for the 
  # denities to plot where x = proportion_success and y = probability
  proportion_success <- c(0, seq(0, 1, length.out = 100), 1)
  dens_curves <- map_dfr(data_indices, function(i) {
    value <- ifelse(i == 0, "Prior", ifelse(data[i], "Success", "Failure"))
    label <- paste0("n=", i)
    probability <- dbeta(proportion_success, 
                         prior_prop[1] + sum(data[seq_len(i)]), 
                         prior_prop[2] + sum(!data[seq_len(i)]))
    probability <- probability / max(probability)
    data_frame(value, label, proportion_success, probability)
  })
  # Turning label and value into factors with the right ordering for the plot
  dens_curves$label <- fct_rev(factor(dens_curves$label, 
  levels =  paste0("n=", data_indices )))
  dens_curves$value <- factor(dens_curves$value, 
  levels = c("Prior", "Success", "Failure"))
  
  p <- ggplot(dens_curves, aes(x = proportion_success, y = label, 
                               height = probability, fill = value)) +
    ggridges::geom_density_ridges(stat="identity", color = "white", alpha = 0.8, 
                                  panel_scaling = TRUE, size = 1) +
    scale_y_discrete("", expand = c(0.01, 0)) +
    scale_x_continuous("Underlying proportion of success") +
    scale_fill_manual(values = hcl(120 * 2:0 + 15, 100, 65), name = "", drop = FALSE,
                      labels =  c("Prior   ", "Success   ", "Failure   ")) +
    ggtitle(paste0(
      "Binomial model - Data: ", sum(data),  " successes, " , 
      sum(!data), " failures")) +
    theme_light() +
    theme(legend.position = "top")
  print(p)  
  
  # Returning a sample from the posterior distribution that can be further 
  # manipulated and inspected
  posterior_sample <- rbeta(n_draws, prior_prop[1] + sum(data), prior_prop[2] + 
    sum(!data))
  invisible(posterior_sample)
}

# simulação dos dados:
big_data <- sample(c(TRUE, FALSE), prob = c(0.75, 0.25), 
                   size = 100, replace = TRUE)

# posterior                   
posterior <- prop_model(big_data)

Warning: `data_frame()` was deprecated in tibble 1.1.0.
ℹ Please use `tibble()` instead.

Warning in ggridges::geom_density_ridges(stat = "identity", color = "white", :
Ignoring unknown parameters: `size`

# intervalo de credibilidade do parâmetro p
quantile(posterior, c(0.025, 0.5, 0.975))

     2.5%       50%     97.5% 
0.6360817 0.7272447 0.8082072 

---

---

---

2. modelos conjugados

As distribuições a priori, ou priores, informativas codificam o conhecimento existente sobre os parâmetros do modelo antes da análise de novos dados.

Em astronomia, as informações a priori frequentemente provêm de levantamentos anteriores, modelos teóricos ou restrições físicas bem estabelecidas.

Vamos agora considerar exemplos que ilustram como diversas distribuições a priori informativas comumente utilizadas surgem naturalmente em problemas astronômicos.

2.1 Prior gaussiano e estimativa da constante de Hubble

O prior gaussiano (ou normal) é adequado quando o valor de um parâmetro é aproximadamente conhecido e queremos considerar novos dados

Exemplo
Considere estimar a constante de Hubble \(H_0\) para um novo conjunto de dados. Estudos anteriores (como o satélite Planck e estudos de distâncias de galáxias) indicam \[ H_0 = 70 \pm 5 ~{\rm km}~{\rm s}^{-1}~{\rm Mpc}^{-1}\] que podemos expressar como um prior \(H_0 \sim N(\mu=70, \sigma = 5)\).

Um novo survey mede \(H_0 = 67 \pm 3~{\rm km}~{\rm s}^{-1}~{\rm Mpc}^{-1}\). Podemos combinar esta nova informação com o prior de \(H_0\) para obter um posterior que considera ambas informações.

Este é um caso de priores conjugados: quando o posterior tem a mesma forma funcional que o prior. Se o prior é \(P(w) \sim N(\mu_0,\sigma_0)\) e a verossimilhança de uma nova medida \(y\) é \(P(y|w) \sim N(w,\sigma)\), o posterior também é gaussiano: \[P(w|y) \sim N(\mu_p,\sigma_p)~~~~~\sigma_p^2 = (1/\sigma_0^2 + 1/\sigma^2)^{-1}~~~~~\mu_p = \sigma_p^2(\mu_0/\sigma_0^2 + y/\sigma^2)\] Logo, a média do posterior é uma média ponderada pela precisão da média do prior e do valor observado. Medidas com menores incertezas têm peso maior.

# bibliotecas gráficas
library(ggplot2)
library(gridExtra)

Attaching package: 'gridExtra'

The following object is masked from 'package:dplyr':

    combine

# vamos precisar de uma função de integração para normalizar o posterior
# vamos usar a regra trapezoidal
trapz <- function(x, y){
  idx <- 2:length(x)
  sum((x[idx] - x[idx-1]) * (y[idx] + y[idx-1]) / 2)
}

# grade de valores de H0 que consideraremos:
# valores entre 45 e 95 km/s/Mpc
H0_grade <- seq(45,95,length=1000)

# prior gaussiano calculado na grade
mu_prior <- 70
sd_prior <- 5

prior <- dnorm(H0_grade, mu_prior, sd_prior)

# nova medida e sua verossimilhança gaussiana: 
H0_obs <- 67
sd_obs <- 3
ver <- dnorm(H0_grade, H0_obs, sd_obs)

# posterior
post <- prior*ver
post <- post / trapz(H0_grade, post)

df1 <- data.frame(
H0 = H0_grade,
prior = prior/max(prior),
verossimilhanca = ver/max(ver),
posterior = post/max(post)
)

df1 <- reshape2::melt(df1,id="H0")

p1 <- ggplot(df1,aes(H0,value,color=variable))+
geom_line(size=1)+
labs(title="prior gaussiano: constante de Hubble",
y="PDF em escala")

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

p1

# estatísticas do posterior:
# moda
H0_grade[which.max(post)]

[1] 67.77277

# largura: desvio padrão
# método: ajustamos uma gaussiana ao log do posterior:
# log(P) = log(A) - (x-mu)^2 / (2*sigma^2)
# isso é uma parábola: y' = a*x^2 + b*x + c, 
# onde a = -1/(2*sigma^2)
# ajuste de um modelo quadrático com lm()
log_y <- log(post)
# remoção de valores -Inf ou NA 
validos <- is.finite(log_y)
model <- lm(log_y[validos] ~ poly( H0_grade[validos], 2, raw=TRUE))

# coeficientes: y = c + b*x + a*x^2
# a é o coeficiente de x^2
a <- coef(model)[3]
# sigma:
sigma_fit <- sqrt(-1 / (2 * a))
sigma_fit

poly(H0_grade[validos], 2, raw = TRUE)2 
                               2.572479 

# valor teórico
sigma_modelo = 1/sqrt(1/sd_prior^2+1/sd_obs^2)
sigma_modelo

[1] 2.572479

---

---

2.2 PRIOR BETA: fração de planetas

A distribuição Beta é amplamente utilizada para parâmetros que representam probabilidades ou frações, pois é definida no intervalo [0, 1]. Vamos designar por \(w\) essas probabilidades ou frações.

Como vimos acima, a distribuição Beta é conjugada da verossimilhança binomial, o que significa que o posterior também será uma distribuição Beta. Se o prior é \(w \sim Beta(a,b)\) e observamos \(k\) sucessos em \(n\) tentativas independentes, a verossimilhança segue uma distribuição binomial \(P(k|n) \sim Binomial(n,w)\) e o posterior será \(P(w|y,n) \sim Beta(a+k,b+n-k)\).

Exemplo:
Suponha que queremos estimar a fração \(w\) de estrelas que hospedam um certo tipo de exoplaneta. Estudos anteriores sugerem que 20% das estrelas possuem exoplanetas desse tipo.

Para uma distribuição Beta(a,b), temos que \(E(w)=a/(a+b)\). Logo, podemos representar o prior como \(w \sim Beta(8,32)\), que tem média \(E[w] = 8/20 = 0.2\).

Um novo levantamento observa 30 estrelas e detecta 9 planetas. A verossimilhança segue uma distribuição binomial, e a distribuição posterior torna-se \(P(w|dados) \sim Beta(17,53)\).

# grade de valores do parâmetro w
w_grade <- seq(0,1,length=1000)

# parâmetros do prior
alpha <- 8
beta <- 32
beta_prior <- dbeta(w_grade,alpha,beta)

# dados observados
n <- 30
k <- 9

# verossimilhança
ver_beta <- dbinom(k,n,w_grade)

# posterior
post_beta <- dbeta(w_grade,alpha+k,beta+n-k)

df2 <- data.frame(
p=w_grade,
prior=beta_prior/max(beta_prior),
posterior=post_beta/max(post_beta)
)

df2 <- reshape2::melt(df2,id="p")

p2 <- ggplot(df2,aes(p,value,color=variable))+
geom_line(size=1)+
labs(title="prior Beta prior: fração de planetas")

p2

---

---

2.3 PRIOR GAMA: taxa de eventos transientes

A distribuição Gama é frequentemente usada para parâmetros que representam taxas positivas, especialmente em processos de Poisson. Vamos representar a taxa por \(\lambda\).

Um prior natural para taxas positivas é a distribuição Gama, \(\lambda \sim Gamma(a,b)\), onde \(a\) e \(b\) são denominados parâmetros de forma e de taxa, respectivamente.

Suponha que observemos \(k\) eventos em um intervalo de tempo \(T\). A verossimilhança é descrita pela distribuição de Poisson, \(P(k|\lambda,T) \sim Poisson(\lambda T)\) e o posterior de \(\lambda\) é uma distribuição Gama: \[P(\lambda|k,T) \sim Gamma(a+k,b+T).\] O parâmetro de forma do posterior aumenta com o número de eventos, enquanto que o parâmetro de taxa aumenta com o tempo de observação. Conforme mais dados são coletados, o posterior fica cada vez mais concentrado. Este modelo é muito usado para estimar taxas de eventos como supernovas, explosões de raios-gama ou detecções de ondas gravitacionais.

Exemplo:
Considere estimar a taxa \(\lambda\) de supernovas em uma dada amostra de galáxias. Estudos anteriores sugerem uma taxa de aproximadamente dois eventos por ano, mas com considerável incerteza. Uma distribuição a priori pode ser escrita como \(P(\lambda) \sim Gamma(4,2)\), que corresponde a uma taxa média de dois eventos por ano.

Suponha que um certo monitoramento observe 12 supernovas em três anos. A verossimilhança é de Poisson e a distribuição a posteriori permanece Gama: \(P(\lambda|dados) \sim Gamma(16,5)\).

A distribuição Gama é definida para quantidades positivas e é conjugada à verossimilhança de Poisson, tornando-a matematicamente conveniente.

# grade de lambdas
lambda <- seq(0.001,20,length=1000)

# prior
shape <- 4
rate <- 2
prior_g <- dgamma(lambda,shape,rate)

# novas observações
# verossimilhança poissoniana
k <- 12
T <- 3
like_g <- dpois(k,lambda*T)

# posterior
post_g <- dgamma(lambda,shape+k,rate+T)

df3 <- data.frame(
lambda=lambda,
prior=prior_g/max(prior_g),
posterior=post_g/max(post_g)
)

df3 <- reshape2::melt(df3,id="lambda")

p3 <- ggplot(df3,aes(lambda,value,color=variable))+
geom_line(size=1)+
labs(title="prior Gama: taxa de transientes")

p3

---

---

---

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

# variáveis:
E = dados[,1]
N = dados[,2]

#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\), a normalização, e \(\beta\), o índice espectral, são os parâmetros do modelo. Isso supõe que os bins são suficientemente estreitos de modo que a taxa dentro de cada um depende apenas de \(E_i\).

Note que todo modelo envolve escolhas!

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

Vamos considerar duas abordagens diferentes.

Método 1: 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
# consideramos intervalos convenientes para a e b
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')

Método 2: esta abordagem é melhor

vamos supor que os bins são independentes e escrever a verossimilhança como

\[ \mathcal{L}(\alpha,\beta) = \prod_{i} \frac{ \left(\alpha E_i^{-\beta}\right)^{N_i} e^{-\alpha E_i^{-\beta}} }{ N_i! }, \] ou \[ \log \mathcal{L}(\alpha,\beta) = \sum_{i} \left[ N_i \log(\alpha E_i^{-\beta}) - \alpha E_i^{-\beta} - \log(N_i!) \right]. \]

Para assegurar a positividade dos parâmetros e melhorar a estabilidade numérica do problema, é comum se fazer uma reparametrização:
\(a = \log \alpha, \qquad b = \log \beta,\) tal que \(\alpha = e^a\) e \(\beta = e^b\).

Vamos adotar priores gaussianos pouco informativos: \(a \sim \mathcal{N}(0, \sigma_a), \qquad b \sim \mathcal{N}(0, \sigma_b),\) com valores altos do desvio padrão (\(\sigma_a = 10\), \(\sigma_b = 5\)).

Nesse caso, o posterior fica: \[ p(a,b \mid \{N_i\}) \propto \mathcal{L}(e^a, e^b) \, \mathcal{N}(a \mid 0, \sigma_a) \, \mathcal{N}(b \mid 0, \sigma_b).\] Vamos amostrar \(a\) e \(b\) por MCMC e depois transformar de volta para \(\alpha\) e \(\beta\):
\(\alpha = e^a, \qquad \beta = e^b\).

# Log-posterior
log_posterior <- function(log_alpha, log_beta, E, N) {

  alpha <- exp(log_alpha)
  beta  <- exp(log_beta)

  mu <- alpha * E^(-beta)

  if (any(!is.finite(mu)) || any(mu <= 0)) {
    return(-Inf)
  }

  # log da verossimilhança poissoniana
  loglik <- sum(dpois(N, lambda = mu, log = TRUE))

  # priores:
  logprior_alpha <- dnorm(log_alpha, mean = 0, sd = 10, log = TRUE)
  logprior_beta  <- dnorm(log_beta,  mean = 0, sd = 5,  log = TRUE)

  loglik + logprior_alpha + logprior_beta
}

# função customizada para MCMC com Metropolis-Hastings 
run_mcmc <- function(E, N,
                     n_iter = 20000,
                     init_alpha = NULL,
                     init_beta = 1,
                     proposal_sd = c(0.08, 0.04)) {

  if (is.null(init_alpha)) {
    init_alpha <- mean(N)
  }

  chain <- matrix(NA, nrow = n_iter, ncol = 2)
  colnames(chain) <- c("alpha", "beta")

  log_alpha <- log(init_alpha)
  log_beta  <- log(init_beta)

  current_lp <- log_posterior(log_alpha, log_beta, E, N)

  accepted <- 0

  for (t in 1:n_iter) {

    prop_log_alpha <- rnorm(1, log_alpha, proposal_sd[1])
    prop_log_beta  <- rnorm(1, log_beta,  proposal_sd[2])

    prop_lp <- log_posterior(prop_log_alpha, prop_log_beta, E, N)

    log_accept_ratio <- prop_lp - current_lp

    if (log(runif(1)) < log_accept_ratio) {
      log_alpha <- prop_log_alpha
      log_beta  <- prop_log_beta
      current_lp <- prop_lp
      accepted <- accepted + 1
    }

    chain[t, ] <- c(exp(log_alpha), exp(log_beta))
  }

  list(
    chain = chain,
    acceptance_rate = accepted / n_iter
  )
}

fit <- run_mcmc(E, N, n_iter = 30000)

chain <- fit$chain
burnin <- 5000
post <- chain[-seq_len(burnin), ]

cat("taxa de aceitação:", fit$acceptance_rate, "\n\n")

taxa de aceitação: 0.3445 

cat("\nsumário do posterior de alpha:\n")

sumário do posterior de alpha:

print(quantile(post[, "alpha"], probs = c(0.025, 0.5, 0.975)))

    2.5%      50%    97.5% 
159.5876 209.7705 279.8110 

cat("\nsumário do posterior de beta:\n")

sumário do posterior de beta:

print(quantile(post[, "beta"], probs = c(0.025, 0.5, 0.975)))

    2.5%      50%    97.5% 
2.373776 2.590292 2.827599 

# plots de diagnóstico
par(mfrow = c(2, 2))

plot(chain[, "alpha"], type = "l",
     xlab = "Iteration", ylab = expression(alpha),
     main = expression("cadeia de " * alpha))

plot(chain[, "beta"], type = "l",
     xlab = "Iteration", ylab = expression(beta),
     main = expression("cadeia de " * beta))

hist(post[, "alpha"], breaks = 40,
     main = expression("Posterior de " * alpha),
     xlab = expression(alpha),col='gold')

hist(post[, "beta"], breaks = 40,
     main = expression("Posterior de " * beta),
     xlab = expression(beta),col='gold')

par(mfrow = c(1, 1))

alpha_hat <- median(post[, "alpha"])
beta_hat  <- median(post[, "beta"])

E_grid <- seq(min(E), max(E), length.out = 200)
N_pred <- alpha_hat * E_grid^(-beta_hat)

plot(E, N, pch = 19,
     xlab = "Energia E",
     ylab = "Countagens N",
     main = "ajuste")

lines(E_grid, N_pred, lwd = 3,col='blue')
lines(E,NC,lwd=3,col='salmon')
legend('topright',legend=c('método 1','método 2'),text.col=c('salmon','blue'),  bty='n')

Compare os resultados dos dois métodos:

# medianas e intervalo de credibilidade de 95%
# método 1:
quantile(p1, probs = c(0.025, 0.5, 0.975))

    2.5%      50%    97.5% 
153.4755 186.0716 223.5588 

quantile(p2, probs = c(0.025, 0.5, 0.975))

    2.5%      50%    97.5% 
2.272523 2.424141 2.567655 

# método 2:
quantile(post[, "alpha"], probs = c(0.025, 0.5, 0.975))

    2.5%      50%    97.5% 
159.5876 209.7705 279.8110 

quantile(post[, "beta"], probs = c(0.025, 0.5, 0.975))

    2.5%      50%    97.5% 
2.373776 2.590292 2.827599 

---

---

---

4. 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 20.80953 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 5.9841 secs

---

---

---

5. 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 1.347018 mins

---

---

---

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?

6.1 abordagem analítica

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

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

6.2 ABC

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

---

---

---

7. Exercício

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

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

Considere os dados

x = c(0.8, 0.4, 1.2, 0.7, 1.8, 0.5, 0.9, 1.4, 0.6, 0.3)

Vamos supor que eles obedeçam a uma distribuição exponencial, \(P(x) \propto \exp(-ax)\), \(a > 0\). Nosso objetivo é determinar \(a\).

Note que o valor médio esperado para \(x\) nesta distribuição é \(1/a\).

1. Qual é o valor de \(a\) obtido da médio dos dados?

2. Vamos usar ABC para estimar \(a\). Para isso:

• 2.1 sorteie os parâmetro \(a\) de uma distribuição uniforme no intervalo \([0, log10(U)]\)

• 2.2 simule um conjunto de dados de mesma dimensão que o observado usando este parâmetro

• 2.3 calcule uma estatística para os dados simulados, como a média, mediana, …

• 2.4 compare esta estatística com a mesma estatística calculada com os dados observados

• 2.5 aceite o parâmetro simulado se as duas estatísticas acima forem suficientemente próximas, menores que um certo \(\epsilon\).

3. Faça 10000 simulações com \(\epsilon\) igual a 0.05 e 0.1. Faça um gráfico do posterior e estime a média do posterior e seu intervalo de credibilidade de 95% para cada caso.

4. Calcule a solução de máxima verossimilhança e compare os resultados.

---

---