#-----------------------------------------------------------------------
# Introdução à inferência Bayesiana 
# Prioris em diferentes parametrizações
#                                            Prof. Dr. Fernando Mayer
#                                       Prof. Dr. Paulo Justiniano R. Jr
#
#                      Laboratory of Statistics and Geoinformation (LEG)
#                Department of Statistics · Federal University of Paraná
#                                        2025-Out-17 · Curitiba/PR/Brazil
#-----------------------------------------------------------------------

## Objetivo: mostrar que definir priori em uma escala (parametrização) implica
## em prioris equivalentes para transformações do parâmetro/reparametrizações

## Ilustrações: 
##  - simular da escala na qual se define a priori
##  - transformar os valores simulados
##  - comparar histograma dos valores simulados transformados com expressões
##    teóricas das densidades na escala transformada obtidas pelo método do
##    Jacobiano de transformação de variáveis aleatórias, quando possível.

##
## Simulação de distribuições a priori
##
par.ori <- par(no.readonly=TRUE)
nsim <- 100000

## Exemplo 1: Definindo priori para proporção $\theta$ ("p") e obtendo as prioris equivalentes
##            em outras escalas/parametrizações
## Y \sim B(\theta)
## \theta \sim B(1, 1) (= U(0,1))
par(mfcol=c(2,3), mar=c(3,3,1,1), mgp=c(2,1,0))
th <- rbeta(nsim, 1, 1)
hist(th, prob=TRUE, 30, xlim = c(0,1), main = ""); lines(density(th))
curve(dbeta(x, 1, 1), from=0, to=1, col=4, add=TRUE)

th2 <- th^2
hist(th2, prob=TRUE, 30, main = ""); lines(density(th2))
## densidade teórica de \theta^2 (por transformação de variáveis)
## J = \frac{1}{2\sqrt{\theta}}
## f(\theta^2) = f(\theta) |J| = 1 * \frac{1}{2\sqrt{\theta^2}} = \frac{1}{2\sqrt{th2}}
##dth2 <- function(th2) 1/(2*sqrt(th2))   # para theta \sim U(0,1)
dth2 <- function(th2, a = 1, b = 1)
    dbeta(sqrt(th2), a, b) * 1/(2*sqrt(th2)) # para theta \sim B(a,b)
curve(dth2, from=0, to=1, col=4, lwd = 2, add=TRUE)

th.odds <- th/(1-th)  ## psi
hist(th.odds, prob=TRUE, 30, main = ""); lines(density(th.odds))
hist(th.odds[th.odds<50], prob=TRUE, 30, main = ""); lines(density(th.odds[th.odds < 50]))
## J = \frac{1}{(1+\psi)^2}
## f(\psi) = f(\theta) |J|
dth.odds <- function(psi, a = 1, b = 1) dbeta(psi/(1+psi), a, b) * 1/((1+psi)^2)
curve(dth.odds, from=0, to=50, col=4, lwd = 2, add=TRUE)

th.logit <- log(th/(1-th))  ## psi = \frac{e^\theta}{1+e^\theta}
hist(th.logit, prob=TRUE, 50, main = ""); lines(density(th.logit))
## J = \frac{\exp(\psi)}{(1+\exp\{\psi\})^2}
## f(\psi) = f(\theta) |J|
dth.logodds <- function(psi, a = 1, b = 1)
    dbeta(exp(psi)/(1+exp(psi)), a, b) * exp(psi)/((1+exp(psi))^2)
curve(dth.logodds, from=-10, to=10, col=4, lwd = 2, add=TRUE)

## Comentário: definir priori em uma parametrização
## implica em prioris "equivalentes" em outras parametrizações

par(mfrow=c(1,3), mar=c(3,3,1,1), mgp=c(2,1,0))
## Mudando parâmetros da Beta: caso "simétrico" com a = b = 5
th <- rbeta(nsim, 5, 5)
hist(th, prob=TRUE, 30, xlim = c(0,1), main = ""); lines(density(th))
curve(dbeta(x, 5, 5), from=0, to=1, col=4, lwd = 2, add=TRUE)

th2 <- th^2
hist(th2, prob=TRUE, 30, main = ""); lines(density(th2))
curve(dth2(x, a = 5, b = 5), from=0, to=1, col=4, lwd = 2, add=TRUE)

th.logit <- log(th/(1-th))
hist(th.logit, prob=TRUE, 30, main = ""); lines(density(th.logit))
curve(dth.logodds(x, 5, 5), from=-10, to=10, col=4, lwd = 2, add=TRUE)

## Beta assimétrica com a = 6, b = 2
par(mfrow=c(1,2), mar=c(3,3,1,1), mgp=c(2,1,0))

th <- rbeta(nsim, 6, 2)
hist(th, prob=TRUE, 30, xlim = c(0,1), main = ""); lines(density(th))
curve(dbeta(x, 6, 2), from=0, to=1, col=4, add=TRUE)

th.logit <- log(th/(1-th))
hist(th.logit, prob=TRUE, 30, main = ""); lines(density(th.logit))
curve(dth.logodds(x, 6, 2), from=-10, to=10, col=4, lwd = 2, add=TRUE)

## priori de Jeffreys: B(1/2, 1/2)
par(mfrow=c(1,3), mar=c(3,3,1,1), mgp=c(2,1,0))

th <- rbeta(nsim, 1/2, 1/2)
hist(th, prob=TRUE, 30, xlim = c(0,1), main = ""); lines(density(th))
curve(dbeta(x, 1/2, 1/2), from=0, to=1, col=4, lwd = 2,  add=TRUE)

th2 <- th^2
hist(th2, prob=TRUE, 30, main = ""); lines(density(th2))
curve(dth2(x, 1/2, 1/2), from=0, to=1, col=4, lwd = 2, add=TRUE)

th.logit <- log(th/(1-th))
hist(th.logit, prob=TRUE, 30, main = ""); lines(density(th.logit))
curve(dth.logodds(x, 1/2, 1/2), from=-10, to=10, col=4, lwd = 2, add=TRUE)

## E se ...? Experimentar outras possibilidades!!!


## Exemplo 2: Definindo prioris para o logito ($\psi$) de uma proporção
## Y \sim B(\theta)
## \psi = logit(\theta) = log(\theta/(1-\theta)) \sim N(mu, \sigma^2)
## \theta = \frac{e^{\psi}}{1+e^{\psi}}

## Distribuição teórica (por transformação de variáveis)
##\frac{d \psi}{d \theta} = \frac{1}{\theta}{1-\theta}
##f_{\theta}(\theta) = f_{\psi}(logito(p); \mu, \sigma) \frac{d \psi}{d \theta} 
##                   = \phi(\log(p/(1-p)); \mu, \sigma) \frac{1}{\theta}{1-\theta} 
##em que \phi(\cdot; \mu, \sigma) é a densidade da distribuição normal de média $\mu$ e desvio padrão $\sigma$.
dinvlogito <- function(x, mean, sd, log = FALSE){
    y <- numeric(length(x))
    IND <- (x > 0 | x < 1)
    x <- x[IND] 
    xlogito <- log(x/(1-x))
    y[IND] <- dnorm(xlogito, mean = mean, sd = sd, log = TRUE) - (log(x)+log(1-x))
    if(!log) y[IND] <- exp(y[IND])
    return(y)
}
par(mfcol=c(1,2), mar=c(3,3,1,1), mgp=c(2,1,0))
curve(dnorm(x, m = 0, sd = 1), from=-4, to=4,
      xlab = expression(psi), ylab = expression(f(psi)))
curve(dinvlogito(x, m=0, sd = 1), from=0, to=1,
      xlab = expression(theta), ylab = expression(f(theta)))
##
## normal para logito
##
par(mfcol=c(2,2), mar=c(3,3,1,1), mgp=c(2,1,0))
th.logit <- rnorm(nsim, m=0, sd=3)
hist(th.logit, prob=TRUE, 30, main = ""); lines(density(th.logit))
curve(dnorm(x, 0, 3), from=-10, to=10, col=4, add=TRUE)

th <- exp(th.logit)/(1+exp(th.logit))
hist(th, prob=TRUE, 30, xlim = c(0,1), col = "transparent", main = ""); lines(density(th))
## densidade teorica
curve(dinvlogito(x, m = 0, sd = 3), from = 0, to = 1, col = 4, lwd = 2, add = TRUE)
## densidade beta aproximada
(bt <- MASS:::fitdistr(th, densfun="beta", start=list(shape1=1, shape2=1))) ## beta que mais se aproxima, mas note que não é a priori correspondente
curve(dbeta(x, coef(bt)[1], coef(bt)[2]), add = TRUE, col = 2)

##
th.logit <- rnorm(nsim, m=0, sd=1)
hist(th.logit, prob=TRUE, 30, main = ""); lines(density(th.logit))
curve(dnorm(x, 0, 1), from=-10, to=10, col=4, add=TRUE)

th <- exp(th.logit)/(1+exp(th.logit))
hist(th, prob=TRUE, 30, xlim = c(0,1), main = ""); lines(density(th))
## densidade teorica
curve(dinvlogito(x, m = 0, sd = 1), from = 0, to = 1, col = 4, lwd = 2, add = TRUE)
## densidade beta aproximada
bt <- MASS:::fitdistr(th, densfun="beta", start=list(shape1=1, shape2=1))
curve(dbeta(x, coef(bt)[1], coef(bt)[2]), add = TRUE, col = 2) # aproxima mas não "bate"

par(mfcol=c(2,4), mar=c(3,3,1,1), mgp=c(2,1,0))
th.logit <- runif(nsim, -100, 100)
hist(th.logit, prob=TRUE, 30, main = ""); lines(density(th.logit))
curve(dunif(x, -100, 100), from=-10, to=10, col=4, add=TRUE)

th <- exp(th.logit)/(1+exp(th.logit))
hist(th, prob=TRUE, 30, xlim = c(0,1), main = ""); lines(density(th))
##
th.logit <- rnorm(nsim, 0, 30)
hist(th.logit, prob=TRUE, 30); lines(density(th.logit))
curve(dnorm(x, 0, 30), from=-10, to=10, col=4, add=TRUE)

th <- exp(th.logit)/(1+exp(th.logit))
hist(th, prob=TRUE, 30, xlim = c(0,1), main = ""); lines(density(th))
## densidade teorica
curve(dinvlogito(x, m = 0, sd = 30), from = 0, to = 1, col = 4, lwd = 2, add = TRUE, n = 500)

##
th.logit <- rnorm(nsim, m=0.5, sd=1)
hist(th.logit, prob=TRUE, 30, main = ""); lines(density(th.logit))
curve(dnorm(x, 0.5, 1), from=-10, to=10, col=4, add=TRUE)

th <- exp(th.logit)/(1+exp(th.logit))
hist(th, prob=TRUE, 30, xlim = c(0,1), main = ""); lines(density(th))
## densidade teórica:
curve(dinvlogito(x, m = 0.5, sd = 1), from = 0, to = 1, col = 4, lwd = 2, add = TRUE)
#bt <- MASS:::fitdistr(th, densfun="beta", start=list(shape1=1, shape2=1))
#curve(dbeta(x, coef(bt)[1], coef(bt)[2]), add = TRUE, col = 2)
##
th.logit <- rnorm(nsim, m=1, sd=2)
hist(th.logit, prob=TRUE, 30, main = ""); lines(density(th.logit))
curve(dnorm(x, 1, 2), from=-10, to=10, col=4, add=TRUE)

th <- exp(th.logit)/(1+exp(th.logit))
hist(th, prob=TRUE, 30, xlim = c(0,1), main = ""); lines(density(th))
## densidade teórica:
curve(dinvlogito(x, m = 1, sd = 2), from = 0, to = 1, col = 4, lwd = 2, add = TRUE)
#bt <- MASS:::fitdistr(th, densfun="beta", start=list(shape1=1, shape2=1))
#curve(dbeta(x, coef(bt)[1], coef(bt)[2]), add = TRUE, col = 2)

##
## Priori Gamma para logito
##
dgammainvlogito <- function(x, shape, rate = 1, scale = 1/rate, log = FALSE){
    y <- numeric(length(x))
    IND <- (x > 0 | x < 1)
    x <- x[IND] 
    xlogito <- log(x)-log(1-x)  # log(x/(1-x))
    y[IND] <- dgamma(xlogito, shape, scale = scale, log = TRUE) - log(x*(1-x))
    if(!log) y[IND] <- exp(y[IND])
    return(y)
}
par(mfrow=c(1,2), mar=c(3,3,1,1), mgp=c(2,1,0))
curve(dgamma(x, shape = 1.5, sc = 1), from=0, to=10, xlab = expression(psi), ylab = expression(f(psi)))
curve(dgammainvlogito(x, shape = 1.5, sc = 0.6), from=0, to=1, xlab = expression(theta), ylab = expression(f(theta)))

par(mfcol=c(2,4), mar=c(3,3,1,1), mgp=c(2,1,0))
th.logit <- rgamma(nsim, sh=1.5, sc=0.6)
hist(th.logit, prob=TRUE, 30, main = ""); lines(density(th.logit))
curve(dgamma(x, sh=1.5, sc=0.6), from=0, to=5, col=4, add=TRUE)

th <- exp(th.logit)/(1+exp(th.logit))
hist(th, prob=TRUE, 30, xlim = c(0,1), main = ""); lines(density(th))
curve(dgammainvlogito(x, shape = 1.5, sc = 0.6), from=0, to=1, col = 4, lwd = 2, add = TRUE)

##
th.logit <- rgamma(nsim, sh=1, sc=1)
hist(th.logit, prob=TRUE, 30, main = ""); lines(density(th.logit))
curve(dgamma(x, sh=1, sc=1), from=0, to=5, col=4, add=TRUE)

th <- exp(th.logit)/(1+exp(th.logit))
hist(th, prob=TRUE, 30, xlim = c(0,1), main = ""); lines(density(th))
curve(dgammainvlogito(x, shape = 1, sc = 1), from=0, to=1, col = 4, lwd = 2, add = TRUE)
#bt <- MASS:::fitdistr(th, densfun="beta", start=list(shape1=1, shape2=1))
#curve(dbeta(x, coef(bt)[1], coef(bt)[2]), add = TRUE, col = 2)
##
th.logit <- rgamma(nsim, sh=2.5, sc=0.6)
hist(th.logit, prob=TRUE, 30, main = ""); lines(density(th.logit))
curve(dgamma(x, sh=2.5, sc=0.6), from=0, to=5, col=4, add=TRUE)

th <- exp(th.logit)/(1+exp(th.logit))
hist(th, prob=TRUE, 30, xlim = c(0,1), main = ""); lines(density(th))
curve(dgammainvlogito(x, shape = 2.5, sc = 0.6), from=0, to=1, col = 4, lwd = 2, add = TRUE)
#bt <- MASS:::fitdistr(th, densfun="beta", start=list(shape1=1, shape2=1))
#curve(dbeta(x, coef(bt)[1], coef(bt)[2]), add = TRUE, col = 2)
##
th.logit <- rgamma(nsim, sh=1.5, sc=1)
hist(th.logit, prob=TRUE, 30, main = ""); lines(density(th.logit))
curve(dgamma(x, sh=1.5, sc=1), from=0, to=5, col=4, add=TRUE)

th <- exp(th.logit)/(1+exp(th.logit))
hist(th, prob=TRUE, 30, xlim = c(0,1), main = ""); lines(density(th))
curve(dgammainvlogito(x, shape = 1.5, sc = 1), from=0, to=1, col = 4, lwd = 2, add = TRUE)
#bt <- MASS:::fitdistr(th, densfun="beta", start=list(shape1=1, shape2=1))
#curve(dbeta(x, coef(bt)[1], coef(bt)[2]), add = TRUE, col = 2)
##

## experimentar outras possibilidades!!!

## Exemplo 3: Priori para parâmetro da Poisson
## Y \sim P(\theta)
## \theta \sim G(shape=\alpha, scale=\beta)
## f(\theta) = \frac{\beta^\alpha}{\Gamma(\alpha)} \theta^{\alpha-1} exp{-\beta\alpha}
par(mfrow=c(1,3), mar=c(3,3,1,1), mgp=c(2,1,0))

th <- rgamma(nsim, 2, 1.5)
hist(th, prob=TRUE, 30, main=expression(theta)); lines(density(th))
curve(dgamma(x, 2, 1.5), from=0, to=8, col=4, add=TRUE)

## Priori equivalente a (probabilidade) de algum evento 
## Tal probabilidade é uma função do parâmetro
## \psi = P(Y > 0) = 1 - P[Y=0] = 1 - exp{-\theta}
psi <- 1 - exp(-th)
hist(psi, prob=TRUE, 30, main=expression(paste(psi ," = P[Y > 0] = 1 - exp(", -theta, ")"))); lines(density(psi))

## \psi = log(theta)
psi <- log(th)
hist(psi, prob=TRUE, 30, main = expression(paste(psi == log(theta)))); lines(density(psi))

## Como simular da priori de Jeffreys neste caso?
## Uma opção é adotar uma priori discreta (aproximar a contínua por discreta)
#par(mfcol=c(1,1), mar=c(3,3,1,1), mgp=c(2,1,0))
#th <- seq(0.1, 5, by=0.1)
#pth <- 1/sqrt(th)
#pth <- pth/sum(pth)
#sum(pth)
#plot(th, pth, type="h")
#sample(th, 10, prob=pth)


## experimentar outras possibilidades!!!

## Exemplo 4: priori na escala log(parâmetro)
## Y \sim P(\theta)
## psi = log(\theta) \sim N(\mu, \sigma^2)

par(mfrow=c(1,2), mar=c(3,3,1,1), mgp=c(2,1,0))
psi <- rnorm(nsim, m=0.3, sd=0.8)
hist(psi, prob=TRUE, 30, main = expression(psi)); lines(density(psi))
curve(dnorm(x, 0.3, 0.8), from=-4, to=4, col=4, add=TRUE)

th <- exp(psi)
hist(th, prob=TRUE, 50, main = expression(theta == exp(psi)), xlim = c(0,30)); lines(density(th))
## teórica neste caso é uma log-normal
curve(dlnorm(x, 0.3, 0.8), from=0, to=8, col=4, add=TRUE)

## Exemplo 5: priori na escala parâmetro prob. de algum evento
## Y \sim P(\theta)
## \theta = log(1-\psi)
## \psi = P(Y > 0) = 1 - P[Y=0] = 1 - exp{-\theta} \sim Beta(3,2)

par(mfrow=c(1,2), mar=c(3,3,1,1), mgp=c(2,1,0))
psi <- rbeta(nsim, 3, 2)
hist(psi, prob=TRUE, 30, main=expression(paste(psi ," = P[Y > 0] = 1 - exp(", -theta, ")"))); lines(density(psi))
curve(dbeta(x, 3, 2), col = 4, add = TRUE)

th <- log(1-psi)
hist(th, prob=TRUE, 30, main=expression(theta)); lines(density(theta))