Estatística Bayesiana

Enunciado dos Exercícios

Exercício 1

De acordo com o enunciado temos os seguintes dados:

• \(p(I) = 0.3\), I = gémeos idênticos
• \(p(GG|I) = 0.5\), GG = duas gémeas raparigas
• \(p(BB|I)\) = 0.5 & BB =dois gémeos rapazes
• \(p(BG|\overline{I})\) = 0.5, BG = gémeos mistos
• \(p(GG|\overline{I})\) = 0.25,
• \(p(BB|\overline{I})\) = 0.25

O enunciado pede \(p(I|GG)\). De acordo com o Teorema de Bayes: \[ p(I|GG)=\frac{p(GG|I) \times p(I)}{p(GG)} \]

Para calcular a evidência de ter duas gémeas: \[ p(GG) = p(GG|I) \times p(I) + p(GG|\overline{I}) \times p(\overline{I}) = 0.5\times 0.3 + 0.25\times 0.7 = 0.325 \]

Assim: \[p(I|GG) = \frac{0.5\times 0.3}{0.325} \approx 0.46\]

Exercício 2

O enunciado dá-nos a seguinte informação:

• \(p(+|G) = 0.983\)
• \(p(-|~G) = 0.988\)
• \(p(G|F) = 0.12\)

F: houve falha no período G: a mulher está grávida +/-: o teste deu positivo/negativo

Assume-se que F e G são independentes, ie, \(F \perp G\). Nas equações seguintes não se refere F, eg, \(p(+) = p(+|F)\)

• 2a) \[p(+) = p(+|G) \times p(G) + p(+|\overline{G}) \times p(\overline{G}) = 0.983 \times 0.12 + (1-0.988) \times (1-0.12) = 0.12852\]

• 2b) \[p(G|+) = \frac{p(+|G) \times p(G)}{p(+)} = \frac{0.983 \times 0.12}{0.12852} = 0.9178338\]

Exercício 3

xs     <- seq(.1,.3,.05)
priors <- c(.25,.35,.2,.15,.05)

prior <- function(x) {
  return (priors[which(xs==x)])
}

likelihood <- function(theta, n.consumers, n.total) {
  return (theta^n.consumers * (1-theta)^(n.total-n.consumers))
}

likelihoods <- likelihood(xs,2,10)

evidence <- sum(priors * likelihoods)

posteriors <- priors * likelihoods / evidence

par(mfrow=c(1,2))
barplot(priors,    ylim=c(0,.6),names.arg=xs,main="Prior")
barplot(posteriors,ylim=c(0,.6),names.arg=xs,main="Posterior")

Em termos numéricos, temos as seguintes mudanças: \[ \begin{array}{ccc} Proporção & à~priori & à~posteriori \\ 0.10 & 0.25 & 0.187 \\ 0.15 & 0.35 & 0.372 \\ 0.20 & 0.2 & 0.233 \\ 0.25 & 0.15 & 0.163 \\ 0.30 & 0.05 & 0.045 \end{array} \]

Se houvessem 20 consumidores em 100 inquiridos, estas seriam as atualizações:

\[ \begin{array}{ccc} Proporção & à~priori & à~posteriori \\ 0.10 & 0.25 & 0.007 \\ 0.15 & 0.35 & 0.335 \\ 0.20 & 0.2 & 0.473 \\ 0.25 & 0.15 & 0.176 \\ 0.30 & 0.05 & 0.009 \end{array} \]

Exercício 4

O enunciado dá-nos a seguinte informação:

• C - Dia de chuva, \(p(C) = 1/3\)

• P - previsão correta de chuva

• \(p(P|C) = 3/4\)

• \(p(P|\overline{C}) = 2/5\)

Queremos saber \(p(C|P)\)

\[p(C|P) = \frac{p(P|C) \times p(C)}{p(P)} = \frac{p(P|C) \times p(C)}{p(P|C) \times p(C) + p(P|\overline{C}) \times p(\overline{C})} = \frac{3/4 \times 1/3}{3/4 \times 1/3 + 2/5 \times 2/3} = 0.4839\]

Exercício 5

linces <- 0:5
priors <- choose(5,linces) * 0.6^linces * 0.4^(5-linces)

likelihood <- function(linces.obs,linces.prior) {
  result = numeric(length(linces.prior))
  for(i in 1:length(linces.prior)) {
    if (linces.prior[i]<linces.obs)
      result[i] <- 0
    else
      result[i] <- (choose(linces.prior[i],linces.obs)
                    *0.3^linces.obs
                    *0.7^(linces.prior[i]-linces.obs))
  }
  return(result)
}

likelihoods <- likelihood(2,linces)

posteriors <- priors * likelihoods / sum(priors * likelihoods)

par(mfrow=c(1,2))
barplot(priors,xlab="número de linces", ylim=c(0,.5),names.arg=linces,  main="Distribuição à Priori")
barplot(posteriors,xlab="número de linces", ylim=c(0,.5),names.arg=linces, main="Dist. após avistar 2 linces")

Exercício 6

• M migração iniciou-se, \(p(M) = 0.4\)
• T tempo em minutos que é preciso esperar

Segundo o enunciado, \(T|M \sim exp(1/20)\) + \(p(T=t|M) = \frac{1}{20} \times e^{-t/20}\) + \(p(T>t|\overline{M}) = 1\), ie, se a migraçao não começou, não haverá avistamentos não interessa o tempo que esperarmos

A CDF da exponencial é \[p(T \lt t|M) = 1-e^{-t/20} \iff p(T>t|M) = e^{-t/20}\]

No R, o CDF da exponencial é calculado pela função pexp: \(p(T \lt t|M)\) é calculado, neste caso, por pexp(t, rate=1/20).

# M = Migração, p(M) = 0.4
p.M <- .4

# p(T>t|M): probabilidade de esperar *T ou mais* minutos (média de 20, ie, lambda = 1/20)
likelihood <- function(t, migracao) {
  if (migracao)
    return (1-pexp(t, rate=1/20))  ## pexp is the cdf of exponential distribution
  return (1)
}

# queremos saber p(~M|T>60). Vamos calcular pelo T. de Bayes:
likelihood(60,FALSE) * (1-p.M) / (likelihood(60,TRUE)*p.M + likelihood(60,FALSE) * (1-p.M))

## [1] 0.9678749

# podemos calcular tb p(M|T>60) e a diferença para um é a nossa resposta
1 - likelihood(60,TRUE) * p.M / (likelihood(60,TRUE)*p.M + likelihood(60,FALSE) * (1-p.M))

## [1] 0.9678749

# Let's plot likelihood(60,TRUE) 

####### Draw the exp distribution

lambda <- 1/20 # sd parameter
lb   <- 60     # lower bound
ub   <- 1000   # upper bound

x  <- seq(0,200,length=1000)
hx <- dexp(x,lambda)
plot(x, hx, type="n", xlab="y", ylab="Density", xlim=c(0,90), ylim=c(0,.035))
i <- x >= lb & x <= ub
lines(x, hx)
# polygon draws the polygons whose vertices are given in x and y
polygon(c(lb,x[i],ub), c(0,hx[i],0), col="red")

area <- pexp(ub, 1/20) - pexp(lb, 1/20)
result <- paste("P( Y >=",lb,") =", signif(area, digits=6))
# place label on top of graph
mtext(result,3)

title(bquote(paste("Exponential Distribution with ", lambda, "= 1/20")))

Exercício 7

xs <- seq(-6,6,.1)
theta1 <- .5
theta2 <- .5
plot(xs, dnorm(xs,1,sqrt(2))*theta1 + dnorm(xs,2,sqrt(2))*theta2, 
     xlab="X",ylab="f.d.p(X)", type="l")

Exercício 8

thetas <- c(1/12, 1/6, 1/4) 
priors <- c(.25,.5,.25)

# esta seria a funçao exacta:

likelihood <- function(n.sixes, theta, throws=1000) {
  result <- 0
  for(i in 1:length(thetas)) {
    result <- result + 
      priors[i] * choose(throws,n.sixes) * 
      thetas[i]^n.sixes * 
      (1-thetas[i])^(throws-n.sixes)
  }
  return (result)
}

# a approx normal da binomial (para n grande, neste caso n=1000)
# é binomial(n,theta) = normal(n.theta, sqrt(n.theta.(1-theta)))

dist.approx <- function(x) {
  n <- 1000
  result <- 0
  for(i in 1:length(thetas)) {
    result <- result + 
      priors[i] * dnorm(x,n*thetas[i],
                        sqrt(n*thetas[i]*(1-thetas[i])))
  }
  return (result)
}
xs <- 1:350
plot(xs,dist.approx(xs),type="l", xlab="num of sixes")
curve(likelihood(x),from=1,to=350,n=350,col="red",add=TRUE)
legend("topright", c("aproximação", "binomial"), col = c("black","red"), lwd=1)

Exercício 9

Alínea 9a)

• \(y\) - número de caras após \(n\) lançamentos
• \(\theta\) probabilidade de sair caras

Temos: + \(\theta \sim Unif(0,1)\) + \(y \sim Binomial(n,\theta)\)

Assim:

\[ \begin{array}{lclr} p(Y=y) & = & \int_0^1 p(Y=y| \theta) p(\theta) d\theta & \\ & = & \int_0^1 p(Y=y| \theta) d\theta & \color{blue}{p(\theta)=1} \\ & = & \int_0^1 \binom n y \theta^y (1-\theta)^{n-y} d\theta & \\ & = & \binom n y \int_0^1 \theta^y (1-\theta)^{n-y} d\theta & \\ & = & \binom n y B(y+1,n-y+1) & \color{blue}{\int_0^1 \frac{1}{B(a,b)}\theta^{a-1}(1-\theta)^{b-1} d\theta = 1} \\ & = & \binom n y \frac{\Gamma(y+1)\Gamma(n-y+1)}{\Gamma(n+2)} & \color{blue}{B(a,b) = \frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}}\\ & = & \frac{n!}{y!(n-y)!} \frac{y!(n-y)!}{(n+1)!} & \color{blue}{\Gamma(n+1) = n!, n \in N}\\ & = & \frac{1}{n+1} & \\ \end{array} \]

Seguem algumas experiências que fiz antes desta resolução que levaram à conclusão empírica que \(p(Y=y|n)=\frac{1}{n+1}\)

# eg:
y <- 3
n <- 6

# analytic solution (found using Mathematica 8)
marginal.dist <- function(y,n) {
  (y<=n)*((choose(n,y) * gamma(n-y+1) * gamma(y+1)) / gamma(n+2))
}
marginal.dist(y,n)

## [1] 0.1428571

# approximation by sum
grid <- 100                 
thetas <- seq(0,1,1/grid)
P <- function(y,n,theta) {  # p(Y=y|theta,n)
  return (choose(n,y) * theta^y * (1-theta)^(n-y))
}
1/grid*sum(P(y,n,thetas))

## [1] 0.1428571

# since all values of theta are equally probable, the event Y=y for diff y's
# is the same, so it only depends of 'n'. Since there are n+1 diff y's to share 
## the probability mass, it should be P(Y=y|n) = 1/(n+1)

marginal.dist2 <- function(y,n) {
  (y<=n)*(1/(n+1))  # avoid if's when using 'outer'
}
  
# plot for all pairs (y,n)
size <- 25; ns <- 1:size; ys <- 1:size
zs <- outer(ys,ns,marginal.dist2)

rotate <- 105;tilt <- 20; parallelness <- 5.0; shadeval <- 0.05; perspcex <- 0.7
persp(ys , ns , zs,
      xlab="y" , ylab="n" , zlab="marginal" , 
      main="" , cex=perspcex, lwd=0.1, 
      xlim=c(0,size) , ylim=c(0,size) , zlim=c(0,1),
      theta=rotate , phi=tilt , d=parallelness, shade=shadeval)
title(bquote(paste("Marginal Distribution P(Y=y|n) = ", frac(1,n+1))))

Alínea 9b)

Agora \(\theta \sim Beta(\alpha,\beta), \alpha,\beta\gt 0\), ie: \[p(\theta) = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} \theta^ {\alpha-1} (1-\theta)^{\beta-1}, \theta \in [0,1]\]

Relembrar que \(y|n \sim Binomial(y;n,\theta)\), ie: \[p(y|\theta) = \binom n y \theta^y (1-\theta)^{n-y}\]

Queremos determinar que existe um \(\lambda, 0 \lt \lambda \lt 1\), tal que: \[E(\theta|y) = \lambda E(\theta) + (1-\lambda)\frac{y}{n}\]

\[ \begin{array}{lclr} p(\theta | y) & = & \frac{p(y|\theta)p(\theta)}{p(y)} & \color{blue}{T.de~Bayes}\\ & = & \frac{\binom n y \theta^y (1-\theta)^{n-y} \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} \theta^ {\alpha-1} (1-\theta)^{\beta-1}}{\int_0^1 \binom n y \theta^y (1-\theta)^{n-y} \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} \theta^ {\alpha-1} (1-\theta)^{\beta-1} d\theta} & \\ & = & \frac{\binom n y \theta^{y+\alpha-1} (1-\theta)^{n-y+\beta-1} }{\frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)} \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} \binom n y \int_0^1 \theta^y (1-\theta)^{n-y} \theta^ {\alpha-1} (1-\theta)^{\beta-1} d\theta} & \\ & = & \frac{\theta^{y+\alpha-1} (1-\theta)^{n-y+\beta-1}}{ \int_0^1 \theta^y (1-\theta)^{n-y} \theta^ {\alpha-1} (1-\theta)^{\beta-1} d\theta} & \\ & = & \frac{ \theta^{y+\alpha-1} (1-\theta)^{n-y+\beta-1}}{ B(y+\alpha,n-y+\beta)} & \end{array} \]

Esta é a PDF de uma \(Beta(y+\alpha,n-y+\beta)\), ou seja, \(\theta|y \sim Beta(y+\alpha,n-y+\beta)\)

O valor esperado é: + \(E(\theta) = \frac{\alpha}{\alpha+\beta}\) + \(E(\theta|y) = \frac{y+\alpha}{n+\alpha+\beta}\)

Temos que encontrar um \(\lambda\) tal que \[\frac{y+\alpha}{n+\alpha+\beta} = \lambda \frac{\alpha}{\alpha+\beta} + (1-\lambda) \frac{y}{n}\]

Após penosas manipulaçoes algébricas chegamos a \(\lambda = \frac{\alpha+\beta}{n+\alpha+\beta}\) que se encontra entre \(]0,1[\).

Alínea 9c)

No caso de que a distribuição à priori de theta seja uniforme, constate que a variância à posteriori de theta é sempre inferior à variância à priori.

\[var(\theta) = 1/12\]

Unif(0,1) é a mesma funçao que Beta(1,1). Assim, pelas regras de atualizaçao: \[\theta|y \sim Beta(y+1,n-y+1)\]

\[ \begin{array}{lclr} var(\theta|y) & = & \frac{(y+1)(n-y+1)}{(\not{y}+1+n-\not{y}+1)^2(\not{y}+1+n-\not{y}+1+1)} & \\ & = & \frac{(y+1)(n-y+1)}{(n+2)^2(n+3)} & \\ & = & \frac{y+1}{n+2} \frac{n-y+1}{n+2} \frac{1}{n+3} & \\ & \leq & \frac{1}{4} \frac{1}{n+3} & \color{blue}{\frac{y+1}{n+2} + \frac{n-y+1}{n+2} = 1, max = 1/2 \times 1/2} \\ & \lt & \frac{1}{4} \frac{1}{3} & \color{blue}{n>1 \Rightarrow n+3\gt 4} \\ & = & \frac{1}{12} & \end{array} \]

Exercício 10

\[ \begin{array}{lclr} p(\theta | x \lt 3) & \propto & p(x \lt 3 | \theta) p(\theta) & \\ & = & [ p(x=0| theta) + p(x=1| theta) + p(x=2| theta) ] p(\theta) & \\ & = & \Big[ \binom {10} 0 \theta^0 + (1-\theta)^{10} + \binom {10} 1 \theta^1 + (1-\theta)^9 + \binom {10} 2 \theta^2 + (1-\theta)^8 \Big] p(\theta) & \color{blue}{x|\theta \sim Binomial(x;10,\theta)} \\ & \propto & \Big[ \binom {10} 0 \theta^0 + (1-\theta)^{10} + \binom {10} 1 \theta^1 + (1-\theta)^9 + \binom {10} 2 \theta^2 + (1-\theta)^8 \Big] (\theta^3 (1-\theta)^3) & \color{blue}{\theta \sim Beta(4,4)} \\ & = & \theta^3(1-\theta)^{13} + 10 \theta^4(1-\theta)^{12} + 45 \theta^5(1-\theta)^{11} & \\ \end{array} \]

xs <- seq(0,1,0.025)
plot(xs,dbeta(xs,4,4), type="l", ylim=c(0,4))

# after n=10 throws there was y<3 heads, ie, y=0,1 or 2 heads
posterior.y.0 <- dbeta(xs,4,14)
posterior.y.1 <- dbeta(xs,5,13)
posterior.y.2 <- dbeta(xs,6,12)

# we give choose(10 y) weight to each posterior
posterior <- (     posterior.y.0 +
              10 * posterior.y.1 +
              45 * posterior.y.2) /56 # the 56 is just to make it smaller

# draw the posterior
points(xs, posterior, type="l", col="red")
legend("topright", c("prior", "posterior"), col = c("black","red"), lwd=1)

xs[which(posterior==max(posterior))]  # where's the mode

## [1] 0.3

Exercício 11

# mean 0.6, sd 0.3
m <- 0.6
s <- 0.3

# the beta parameters are:
alfa <- m*(m*(1-m)/s^2 - 1)
beta <- (1-m)*(m*(1-m)/s^2 - 1)

xs <- seq(0,1,0.0005)
plot(xs,dbeta(xs,alfa,beta), type="l", xlim=c(.5,1), ylim=c(0,30))

n <- 1000
y <- 650

posterior <- dbeta(xs,alfa+y,beta+(n-y))
points(xs, posterior, type="l", col="red")
legend("topleft", c("prior", "posterior"), col = c("black","red"), lwd=1) 

post.mean <- (y+alfa)/(n+alfa+beta)
post.var <- ((alfa+y)*(beta+n-y))/((alfa+beta+n)^2*(alfa+beta+n+1))

text(.9, 21, paste("mean =",signif(m,3)))
text(.9, 19, paste("var =",signif(s,3)))
text(.9, 17, paste("mean|y =",signif(post.mean,3)))
text(.9, 15, paste("var|y =",signif(post.var,3)))

Exercício 12

É dada uma amostra de \(x=34, n=70\). Considerou-se que a distribuição à priori em relaçao à probabilidade \(\theta\) de sobreviver à doença tem média \(0.4\) e variância \(0.02\).

x <- 34
n <- 70
pre.mean <- .4
pre.var  <- .02
pre.sd   <- sqrt(pre.var)
pre.sd

## [1] 0.1414214

get.parameters.beta <- function(beta.mean, beta.sd) {
  a <- beta.mean*(beta.mean*(1-beta.mean)/beta.sd^2 - 1)
  b <- (1-beta.mean)*(beta.mean*(1-beta.mean)/beta.sd^2 - 1)
  c(a,b)
}

# finding which beta parameters have this statistics
pre.params <- get.parameters.beta(pre.mean,pre.sd)
pre.alfa <- pre.params[1]  
pre.beta <- pre.params[2]  

# plot the beta
thetas <- seq(0, 1, 0.01)
plot(thetas, dbeta(thetas, pre.alfa, pre.beta), type = "l", xlab="theta", ylab="density", main="Prior")
mtext(paste("alfa =",pre.alfa,", beta = ", pre.beta),3)

Qual a probabilidade à priori que 50% dos individuos sobrevivam à doença?

1-pbeta(.5,pre.alfa,pre.beta)

## [1] 0.2445096

Como se afirma que a distribuição à priori é uma \(Beta(4.4, 6.6)\), e a amostra pode ser modelada por uma binomial, teremos também uma distribuição à posteriori \(Beta(\alpha^{*},\beta^{*}) = Beta(\alpha + x, \beta + (n-x))\).

pos.alfa <- pre.alfa + x
pos.beta <- pre.beta + (n-x)
pos.mean <- pos.alfa / (pos.alfa+pos.beta)
pos.mean # the mean

## [1] 0.4740741

pos.var  <- (pos.alfa*pos.beta) / ((pos.alfa+pos.beta)^2*(pos.alfa+pos.beta+1))
pos.var # the variance

## [1] 0.003040583

pos.sd   <- sqrt(pos.var)
pos.sd # the st.deviation

## [1] 0.05514149

pos.mode <- (pos.alfa-1) / (pos.alfa+pos.beta-2)
pos.mode # the mode

## [1] 0.4734177

Qual a probabilidade à posteriori que 50% dos individuos sobrevivam à doença?

1-pbeta(.5,pos.alfa,pos.beta)

## [1] 0.319469

# plot the beta
plot(thetas, dbeta(thetas, pos.alfa, pos.beta), type = "l", xlab="theta", ylab="density", main="Prior vs Posterior")
points(thetas, dbeta(thetas, pre.alfa, pre.beta), type = "l", col="red")
mtext(paste("pos.alfa =",pos.alfa,", pos.beta = ", pos.beta),3)
legend("topright", c("prior", "posterior"), col = c("red","black"), lwd=2)

São os dados informativos em relaçao a \(theta\)?

Sim, há uma grande diferença entre as duas distribuições, antes e depois de considerar os dados.

Assuma agora que se trata de uma doença pouco frequente e que o médico não consegue exprimir a sua opinião acerca de \(theta\), antes de observar os dados. Qual será então a distribuição à priori para \(theta\) que deverá considerar? Que modificações se observam na distribuição à posteriori?

Dever-se-á usar uma beta uniforme:

pre.alfa.2 <- 1
pre.beta.2 <- 1
plot(thetas, dbeta(thetas, pre.alfa.2, pre.beta.2), type = "l", xlab="theta", ylab="density", main="Prior")
mtext(paste("alfa =",pre.alfa.2,", beta = ", pre.beta.2),3)

Se fizermos os mesmos cálculos com os mesmos dados:

pos.alfa <- pre.alfa.2 + x
pos.beta <- pre.beta.2 + (n-x)
pos.mean <- pos.alfa / (pos.alfa+pos.beta)
pos.mean # the mean

## [1] 0.4861111

pos.var  <- (pos.alfa*pos.beta) / ((pos.alfa+pos.beta)^2*(pos.alfa+pos.beta+1))
pos.var # the variance

## [1] 0.003422015

pos.sd   <- sqrt(pos.var)
pos.sd # the st.deviation

## [1] 0.05849799

pos.mode <- (pos.alfa-1) / (pos.alfa+pos.beta-2)
pos.mode # the mode

## [1] 0.4857143

Qual a probabilidade à posteriori que 50% dos individuos sobrevivam à doença?

1-pbeta(.5,pre.alfa.2,pre.beta.2)

## [1] 0.5

1-pbeta(.5,pos.alfa,pos.beta)

## [1] 0.4062943

# plot the beta
plot(thetas, dbeta(thetas, pos.alfa, pos.beta), type = "l", xlab="theta", ylab="density", main="Prior vs Posterior")
points(thetas, dbeta(thetas, pre.alfa.2, pre.beta.2), type = "l", col="red")
mtext(paste("pos.alfa =",pos.alfa,", pos.beta = ", pos.beta),3)
legend("topright", c("prior", "posterior"), col = c("red","black"), lwd=2)

A resposta em relação à probabilidade de sobrevivência de 50% dos pacientes mudou menos que no caso anterior.

Exercício 13

O número de defeitos na fita segue uma distribuição de Poisson(\(theta\)), \(theta>0\). A distribuiçao à priori de \(theta\) é uma Gama(3,1).

Em cinco rolos de fita magnética observaram-se os seguintes números de defeitos: 2, 2, 6, 0 e 3, respectivamente para os rolos 1, 2, …,5.

pre.a <- 3
pre.b <- 1

sum.xs <- sum(c(2,2,6,0,3))
n      <- 5
pos.a  <- pre.a + sum.xs
pos.b  <- pre.b + n

thetas <- seq(1,10,.1)
plot(thetas, dgamma(thetas, shape=pre.a, rate=pre.b), type = "l", xlab="theta", ylab="density", main="Prior vs Posterior", col="red", ylim=c(0,2/3))
points(thetas, dgamma(thetas, shape=pos.a, rate=pos.b), type = "l", col="black")
legend("topright", c("prior", "posterior"), col = c("red","black"), lwd=2)
mtext(paste("pos.a =",pos.a,", pos.b = ", pos.b),3)

Exercício 14

1. Supondo que as distribuições à priori de \(\theta_1\) e \(\theta_2\) são \(Gama(2,1)\), obtenha as distribuições à posteriori de \(\theta_1\) e de \(\theta_2\).

2. Represente graficamente as distribuições à posteriori de \(\theta_1\) e de \(\theta_2\) obtidas na alínea anterior. Indique os valores da média e da moda à posteriori.

3. Obtenha os intervalos de credibilidade a 95%, com iguais probabilidades de cauda, para \(\theta_1\) e \(\theta_2\) à posteriori.

pre.a <- 2  # trying not to be informative: has mean 2 and variance 2
pre.b <- 1
thetas <- seq(0,4,.01)
plot(thetas, dgamma(thetas, shape=pre.a, rate=pre.b), type = "l", xlab="theta", ylab="density", col="black", ylim=c(0,3))

# for non educated women
pos.a.ned <- pre.a + 217
pos.b.ned <- pre.b + 111

# for educated women
pos.a.ed <- pre.a + 66
pos.b.ed <- pre.b + 44

# ploting all together
points(thetas, dgamma(thetas, shape=pos.a.ned, rate=pos.b.ned), type = "l", col="red")
points(thetas, dgamma(thetas, shape=pos.a.ed, rate=pos.b.ed), type = "l", col="blue")
legend("topright", c("prior", "post non educated", "post educated"), col = c("black","red", "blue"), lwd=2)

# mode for theta_1 (non educated)
(pos.a.ned-1) / (pos.b.ned)

## [1] 1.946429

The typical statistics for a \(Gamma(\alpha,\beta)\) distribution are:

• Mode = \(\frac{\alpha-1}{\beta}\)
• Mean = \(\frac{\alpha}{\beta}\)

# mode and mean for theta_1 (non educated)
(pos.a.ned-1) / pos.b.ned

## [1] 1.946429

pos.a.ned / pos.b.ned

## [1] 1.955357

# mode and mean for theta_2 (educated)
(pos.a.ed-1) / pos.b.ed

## [1] 1.488889

pos.a.ed / pos.b.ed

## [1] 1.511111

Resta a alínea c) Obtenha os intervalos de credibilidade a 95%, com iguais probabilidades de cauda, para \(\theta_1\) e \(\theta_2\) à posteriori.

Podemos calcular estes intervalos encontrando o percentil 2.5 e 97.5 das respectivas distribuições. Para tal usamos a função qgamma

# 95% confidence interval for non educated women
qgamma(c(.025,.975),shape=pos.a.ned, rate=pos.b.ned)

## [1] 1.704943 2.222679

# 95% confidence interval for educated women
qgamma(c(.025,.975),shape=pos.a.ed, rate=pos.b.ed)

## [1] 1.173437 1.890836

De acordo com esta página o menor intervalo de credibilidade, designado por HPD (Highest Posterior Density) para uma Gama pode ser aproximado por uma normal com média \(\overline{x}\) e variância \(\frac{\overline{x}}{n}\), ie,

\[\Big[ \Big( \overline{x} - 1.96 \sqrt{\frac{\overline{x}}{n}} \Big), \Big( \overline{x} + 1.96 \sqrt{\frac{\overline{x}}{n}} \Big) \Big]\]

mean.ned <- 1.95
n.ned <- 111
paste("non educated 95% HPD: [", signif(mean.ned - 1.96*sqrt(mean.ned/n.ned),5),
      ",", signif(mean.ned + 1.96*sqrt(mean.ned/n.ned),5),"]")

## [1] "non educated 95% HPD: [ 1.6902 , 2.2098 ]"

mean.ed <- 1.5
n.ed <- 44
paste("educated 95% HPD: [", signif(mean.ed - 1.96*sqrt(mean.ed/n.ed),5),
      ",", signif(mean.ed + 1.96*sqrt(mean.ed/n.ed),5),"]")

## [1] "educated 95% HPD: [ 1.1381 , 1.8619 ]"

Exercício 15

Tempo de vida dado por \(f(x|\theta) = \theta e^{-\theta x}, x \gt 0\), sendo \(\theta \sim Gama(\alpha,\beta)\)

Coeficiente de variância = \(\frac{\sigma}{\mu} = 0.5\)

O \(\sigma\) da Gama é \(\frac{\sqrt{\alpha}}{\beta}\) e \(\mu\) é \(\frac{\alpha}{\beta}\), ou seja,

\[\frac{\frac{\sqrt{\alpha}}{\beta}}{\frac{\alpha}{\beta}} = \frac{\sqrt{\alpha}}{\alpha} = \frac{1}{\sqrt{\alpha}} = 0.5 \iff \alpha = 4\]

A posterior tem de ser tal que o CV seja \(0.1\), ou seja \(\alpha = 100\)

Pela regra da atualização da posterior, o acréscimo de \(100-4=96\) no hiperparâmetro \(\alpha\) é a contribuição da dimensão da amostra. Assim \(n = 96\).

Exercício 16

\(f(x|\theta) = \theta e^{-\theta x}, x,\theta \gt 0\)

\(n=20, \overline{x} = 3.8\)

Uma distribuiçao gama de média 0.2 e desvio padrão 0.1 corresponde a \(\theta \sim Gama(\frac{1}{25}, \frac{1}{5})\)

pre.a <- 1/25
pre.b <- 1/5
thetas <- seq(0,4,.01)
plot(thetas, dgamma(thetas, shape=pre.a, rate=pre.b), type = "l", xlab="theta", ylab="density", col="black", ylim=c(0,3))

pos.a <- pre.a + 3.8*20
pos.b <- pre.b + 20
points(thetas, dgamma(thetas, shape=pos.a.ed, rate=pos.b.ed), type = "l", col="red")
legend("topright", c("prior", "post"), col = c("black","red"), lwd=2)

text(3.5, 2, paste("mean =",signif(pre.a/pre.b,3)))
text(3.5, 1.75, paste("var =",signif(pre.a/pre.b^2,3)))
text(3.5, 1.5, paste("mean|x =",signif(pos.a/pos.b,3)))
text(3.5, 1.25, paste("var|x =",signif(pos.a/pos.b^2,3)))

Exercício 17

Sabe-se que \[\frac{1}{Var(\mu|x)} = \frac{n}{\sigma^2} + \frac{1}{d^2}\] onde \(\sigma^2 = 72\) é a variância dos dados e \(d^2=12\) é a variância da prior. Assim, \[\frac{1}{Var(\mu|x)} = \frac{n}{72} + \frac{1}{12} \iff Var(\mu|x) = \frac{72}{n+6}\]

Para que \(Var(\mu|x) \leq 1\), o valor de \(n\) não pode ser inferior a \(66\).

Exercício 18

\(\theta \sim N(b,10^2)\)

\(\theta|x \sim N(52,\sqrt{10}^2)\)

\(x = (x_1,\ldots,x_4)\) com \(\overline{x}=55\)

Sabe-se que \(\theta|x \sim N\Big( \frac{n\tau\overline{x}+cb}{c+n\tau}, \frac{1}{c+n\tau} \Big) = N(52,10)\) onde \(n=4, \tau=\frac{1}{\sigma^2}=\frac{1}{100}\). Temos assim:

\[\left\{ \begin{array}{l} \frac{220/100+bc}{c+4\frac{1}{100}} = 52 \\ \frac{1}{c+4\frac{1}{100}} = 10 \end{array} \right. \]

Resolvendo obtemos \(b=50, c=\frac{3}{50}\) logo \(d^2=\frac{1}{c}\iff d=4.08\).

pre.m  <- 50
pre.sd <- 10
thetas <- seq(-10,80,.01)
plot(thetas, dnorm(thetas, pre.m, pre.sd), type = "l", xlab="theta", ylab="density", col="black", ylim=c(0,.15))

pos.m <- 52
pos.sd <- sqrt(10)
points(thetas, dnorm(thetas, pos.m, pos.sd), type = "l", col="red")

data.m <- 55
data.sd <- 4.08
points(thetas, dnorm(thetas, data.m, data.sd), type = "l", col="blue")
legend("topright", c("prior", "post", "data"), col = c("black","red", "blue"), lwd=2)

Exercício 19

\(\theta \sim N(a,b^2)\)

\(a = E(\theta) = 6\)

Sabemos que a distribuição normal \(N(a,b^2)\) (cuja variância ainda não sabemos, já \(\mu=a=6\)) terá de satisfazer \(P(\theta \gt 12) = 0.15 \iff P(\theta \leq 12) = 0.85\).

Para a distribuição normal padrão \(N(0,1)\):

• pdf: \(\phi(x) = \frac{1}{\sqrt{2\pi}} e^{-\frac{x^2}{2}}\)
• cdf: \(\Phi(x) = \int_{-\infty}^x \frac{1}{\sqrt{2\pi}} e^{-\frac{t^2}{2}} dt\)

Para uma distribuição normal \(N(\mu,\sigma^2)\):

• pdf: \(f(x) = \frac{1}{\sigma} \phi(\frac{\mu-x}{\sigma})\)
• cdf: \(F(x) = \Phi(\frac{x-\mu}{\sigma})\). Esta função é dada por pnorm()
• inversa cdf: \(F^{-1}(p) = \mu + \sigma \Phi^{-1}(p)\), com \(p \in [0,1]\). Esta função é dada por qnorm()

Voltando ao exercício, dizer que \(P(\theta \leq 12) = 0.85\) é o mesmo que \(F(12) = 0.85\). Ora,

\[F(12) = 0.85 \iff \Phi\Big(\frac{12-a}{b}\Big) = \Phi\Big(\frac{12-6}{b}\Big) = \Phi\Big(\frac{6}{b}\Big) = 0.85\]

\[\frac{6}{b} = \Phi^{-1}(0.85) \approx 1.036 \iff b = 5.79\]

Assim, \(\theta \sim N(6,5.79^2)\)

Exercício 20

\(\theta \sim N(80,18^2)\)

\(X \sim N(\theta,10^2)\)

# n: number of obs for X
# x.mean: average of obs
# sigma2: variance of X
# b: mean of prior
# d2: variance of prior
# returns a vector with (mean of posterior, variance of posterior)
find.posterior.param.normal <- function(n, x.mean, sigma2, b, d2) {
  tau <- 1/sigma2
  c <- 1/d2
  new.mu     <- (n*tau*x.mean + c*b) / (c + n*tau) 
  new.sigma2 <- 1/(c + n*tau) 
  return (c(new.mu, new.sigma2))
}

find.posterior.param.normal(50,68,10^2,80,18^2)

## [1] 68.07362  1.98773

So, the posterior distribution is \(N(68,1.41^2)\)

A distribuição preditiva dada a posterior é:

\[y^* \sim N(68, 1.41^2 + 10^2) = N(68,10.1^2)\]

Exercício 21

Qual a distribuição normal à priori sabendo que 25% não falam aos 8 meses e que apenas 9% excedem os 20 meses sem falar? Indique a posterior sabendo que os dados provêm de uma normal com desvio padrão \(\sqrt{63}\)

Para determinar os parâmetros da normal usamos a função get.norm.par() do pacote rriskDistributions. Uma outra função interessante é a fitdistr() do pacote MASS mas que encontra os parametros dada uma amostra.

library(rriskDistributions)

prior.dist <- get.norm.par(p=c(.25,.91),q=c(8,20), plot=FALSE)

## The fitting procedure 'L-BFGS-B' was successful! 
## $par
## [1] 12.016265  5.954611
## 
## $value
## [1] 3.16556e-12
## 
## $counts
## function gradient 
##        7        7 
## 
## $convergence
## [1] 0
## 
## $message
## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

prior.m <- prior.dist["mean"]
prior.sd <- prior.dist["sd"]

pnorm(8,prior.m,prior.sd)

## [1] 0.2500032

1-pnorm(20,prior.m,prior.sd)

## [1] 0.08999835

thetas <- seq(0,30,0.01)
plot(thetas,dnorm(thetas,prior.m,prior.sd),type="l",ylim=c(0,.25))

data <- c(15,26,10,9,15,20,18,11,8,20,7,9,10,11,11,10,12,24,17,11,10)

posterior.params <- find.posterior.param.normal(length(data),mean(data),63,prior.m,prior.sd^2)
points(thetas,dnorm(thetas,posterior.params[1],sqrt(posterior.params[2])),type="l",col="red")
legend("topright", c("prior", "post"), col = c("black","red"), lwd=2)

Qual a probabilidade de \(\theta\) ser superior a 24 meses?

# prior
1 - pnorm(24,prior.m,prior.sd)

## [1] 0.02208291

# posterior
1 - pnorm(24,posterior.params[1],sqrt(posterior.params[2]))

## [1] 9.462364e-11

E se tivéssemos usado uma prior vaga?

prior.m = 0
prior.sd = 1e4 # valores tipicos para uma normal vaga
plot(thetas,dnorm(thetas,prior.m,prior.sd),type="l",ylim=c(0,.25))
posterior.params <- find.posterior.param.normal(length(data),mean(data),63,prior.m,prior.sd^2)
points(thetas,dnorm(thetas,posterior.params[1],sqrt(posterior.params[2])),type="l",col="red")
legend("topright", c("prior", "post"), col = c("black","red"), lwd=2)

# prior
1 - pnorm(24,prior.m,prior.sd)

## [1] 0.4990425

# posterior
1 - pnorm(24,posterior.params[1],sqrt(posterior.params[2]))

## [1] 7.313129e-10

Não faria muita diferença…

Se observarmos mais 10 crianças, com as condições da primeira prior, qual a distribuição preditiva para a idade média de quando irão falar pela primeira vez?

prior.dist <- get.norm.par(p=c(.25,.91),q=c(8,20), plot=FALSE)

## The fitting procedure 'L-BFGS-B' was successful! 
## $par
## [1] 12.016265  5.954611
## 
## $value
## [1] 3.16556e-12
## 
## $counts
## function gradient 
##        7        7 
## 
## $convergence
## [1] 0
## 
## $message
## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

prior.m <- prior.dist["mean"]
prior.sd <- prior.dist["sd"]

posterior.params <- find.posterior.param.normal(length(data),mean(data),63,prior.m,prior.sd^2)

# predictive mean
posterior.params[1]

##       sd 
## 13.40621

# predictive sd
sqrt(posterior.params[2]^2+63/10)

##       sd 
## 3.735052

plot(thetas,dnorm(thetas,posterior.params[1], sqrt(posterior.params[2]^2+63/10)),type="l")

Exercício 22

Seja \(X_1, \ldots X_n\) um conjunto de \(n\) v.a. iid com distribuição uniforme no intervalo \([ 0,\theta ]\) \[f(x_i|\theta) = \frac{1}{\theta}, 0 \lt x_i \lt \theta\]

Assuma que a distribuição à priori para \(\theta\) é a seguinte distribuição imprópria \[p(\theta) = \frac{1}{\theta}, \theta \gt 0\]

1. Prove que a distribuição à posteriori tem a seguinte expressão \[p(\theta|x) = \frac{nM^n}{\theta^{n+1}}, \theta \geq M\] sendo \(M = \operatorname*{arg\,max}_i ~ x_i\).

Resposta: Sabemos que a posterior é dada por: \[p(\theta|x) = \frac{f(x|\theta) p(\theta)}{p(x)}\]

Como as variáveis \(X_i\) são independentes entre si: \[f(x|\theta) = \prod_i f(x_i|\theta) = \frac{1}{\theta^n}\]

Podemos agora calcular a evidência: \[p(x) = \int_\Theta f(x|\theta) p(\theta) ~ d\theta\]

Como para cada \(x_i\) é necessário satisfazer \(0 \lt x_i \lt \theta\), os valores admissíveis do \(\theta\) têm de ser maiores que todos os valores de \(x_i\), i.e., \(\theta \gt M\). Assim:

\[ \begin{array}{lclr} p(x) & = & \int_M^{+\infty} f(x|\theta) p(\theta) ~ d\theta & \\ & = & \int_M^{+\infty} \frac{1}{\theta^n} \frac{1}{\theta} ~ d\theta & \color{blue}{p(\theta)=\frac{1}{\theta}} \\ & = & \int_M^{+\infty} \frac{1}{\theta^{n+1}} ~ d\theta & \\ & = & \Big[ -\frac{1}{n} \frac{1}{\theta^n} \Big]_M^{+\infty} & \\ & = & -\frac{1}{n} \Big( \lim\nolimits_{\theta\to \infty} \tfrac{1}{\theta^n} - \frac{1}{M^n} \Big) & \\ & = & -\frac{1}{n} ( 0 - \frac{1}{M^n} ) & \\ & = & \frac{1}{n \times M^n} & \\ \end{array} \]

Temos calculadas todas as componentes, podemos encontrar a expressão da posterior: \[p(\theta|x) = \frac{1/\theta^{n+1}}{1/nM^n} = \frac{nM^n}{\theta^{n+1}}\]

2. Prove que \(E(\theta|x) = \frac{n}{n-1}M\)

\[ \begin{array}{lclr} E(\theta|x) & = & \int_\Theta \theta p(\theta|x) ~ d\theta & \\ & = & \int_M^{+\infty} \theta \frac{nM^n}{\theta^{n+1}} ~ d\theta & \\ & = & n M^n \int_M^{+\infty} \frac{1}{\theta^n} ~ d\theta & \\ & = & n M^n \frac{1}{n-1} \frac{1}{M^{n-1}} & \\ & = & \frac{n}{n-1}M & \\ \end{array} \]

3. Calcule o valor de \(p(\theta \gt tM | x)\) para \(t \geq 1\)

\(p(\theta \gt tM | x) = 1 - p(\theta \leq tM | x)\)

Sabemos a pdf \[f(\theta|x) = \frac{nM^n}{\theta^{n+1}}\]

Pretendemos calcular a cdf

\[ \begin{array}{lclr} F(\theta|x) & = & \int_M^{\theta} f(\theta|x) ~ d\theta & \\ & = & \int_M^{\theta} \frac{nM^n}{\theta^{n+1}} ~ d\theta & \\ & = & M^n \times (M^{-n} - \theta^{-n}) & \\ & = & 1 - \frac{M^n}{\theta^n} & \\ \end{array} \]

Assim: \[p(\theta \gt tM | x) = 1 - p(\theta \leq tM | x) = 1 - 1 + \frac{M^n}{(tM)^n} = \frac{M^n}{(tM)^n} = \frac{1}{t^n}\]

Exercício 23

Obtenha a distribuição à priori de Jeffreys para:

1. a Poisson(\(\theta\)), \(\theta \gt 0\)

A distribuição de Poisson é dada por \[p(x|\theta) = \frac{\theta^x}{x!} e^{-\theta}\]

O logaritmo da verosimilhança é: \[log(p(x|\theta)) = L = x log(\theta) - \theta - log~x!\]

A sua primeira derivada: \[\frac{\delta L}{\delta\theta} = \frac{x-\theta}{\theta}\]

E a segunda derivada: \[\frac{\delta^2 L}{\delta\theta^2} = -\frac{1}{\theta} + \frac{x-\theta}{\theta^2}\]

Assim:

\[ \begin{array}{lclr} I(\theta) & = & - E_{X|\theta} \Big[ \frac{\delta^2 L}{\delta\theta^2} \Big] & \\ & = & - E_{X|\theta} \Big[ -\frac{1}{\theta} + \frac{X-\theta}{\theta^2} \Big] & \\ & = & - \Big[ -\frac{1}{\theta} + \frac{\theta-\theta}{\theta^2} \Big] & \color{blue}{E_{X|\theta}(X) = \theta} \\ & = & - \Big[ -\frac{1}{\theta} \Big] & \\ & = & \frac{1}{\theta} & \\ \end{array} \]

A prior de Jeffreys é dada pela raiz quadrada deste valor: \[p(\theta) = \sqrt{I(\theta)} = \sqrt{\frac{1}{\theta}}\]

Esta prior é imprópria (não é uma densidade, o integral dá infinito). Seria uma \(gamma(1/2,0)\) se o 2º parametro não tivesse de ser positivo. Porém, isto dá-nos uma indicação para a escolha de uma distribuição, eg, \(gamma(1/2,0.0001)\) que se aproxima da prior de Jeffreys mas é própria!

---

2. a normal N(\(\mu,\sigma^2\)) com \(\mu \in R\) e \(\sigma \gt 0\) conhecido.

A distribuição da normal é dada por \[p(x|\mu) = \frac{1}{\sqrt{2 \pi \sigma^2}} e^{-\frac{(x-\mu)^2}{2\sigma^2}}\]

O logaritmo da verosimilhança é: \[log(p(x|\mu)) = L = log \frac{1}{\sqrt{2 \pi \sigma^2}} - \frac{(x-\mu)^2}{2\sigma^2} \]

A sua primeira derivada: \[\frac{\delta L}{\delta\mu} = \frac{(x-\mu)}{\sigma^2}\]

E a segunda derivada: \[\frac{\delta^2 L}{\delta\mu^2} = - \frac{1}{\sigma^2}\]

Assim:

\[ \begin{array}{lclr} I(\mu) & = & - E_{X|\mu} \Big[ \frac{\delta^2 L}{\delta\mu^2} \Big] & \\ & = & - E_{X|\mu} \Big[ - \frac{1}{\sigma^2} \Big] & \\ & = & - \Big[ - \frac{1}{\sigma^2} \Big] & \color{blue}{E_{X|\sigma}(a) = a} \\ & = & \frac{1}{\sigma^2} & \\ \end{array} \]

A prior de Jeffreys é dada pela raiz quadrada deste valor: \[p(\mu) = \sqrt{I(\mu)} = \sqrt{\frac{1}{\sigma^2}} = \frac{1}{\sigma} \propto 1\]

ou seja, é a uniforme, dado que \(\sigma\) é fixo (é uma constante dada).

---

3. a normal N(\(\mu,\sigma^2\)) com \(\mu \in R\) conhecido e \(\sigma \gt 0\).

A distribuição da normal é dada por \[p(x|\sigma) = \frac{1}{\sqrt{2 \pi \sigma^2}} e^{-\frac{(x-\mu)^2}{2\sigma^2}}\]

O logaritmo da verosimilhança é: \[log(p(x|\sigma)) = L = log \frac{1}{\sqrt{2 \pi \sigma^2}} - \frac{(x-\mu)^2}{2\sigma^2} \]

A sua primeira derivada: \[\frac{\delta L}{\delta\sigma} = - \frac{1}{\sigma} + \frac{(x-\mu)^2}{\sigma^3}\]

E a segunda derivada: \[\frac{\delta^2 L}{\delta\sigma^2} = \frac{1}{\sigma^2} - \frac{3(x-\mu)^2}{\sigma^4}\]

Assim:

\[ \begin{array}{lclr} I(\sigma) & = & - E_{X|\sigma} \Big[ \frac{\delta^2 L}{\delta\sigma^2} \Big] & \\ & = & - E_{X|\sigma} \Big[ \frac{1}{\sigma^2} - \frac{3(X-\mu)^2}{\sigma^4} \Big] & \\ & = & - \Big[ \frac{1}{\sigma^2} - \frac{3\sigma^2}{\sigma^4} \Big] & \color{blue}{E_{X|\sigma}((X-\mu)^2 = \sigma^2} \\ & = & \frac{2}{\sigma^2} & \\ & \propto & \frac{1}{\sigma^2} & \\ \end{array} \]

A prior de Jeffreys é dada pela raiz quadrada deste valor: \[p(\sigma) = \sqrt{I(\sigma)} = \sqrt{\frac{1}{\sigma^2}} = \frac{1}{\sigma}\]

---

4. a geométrica dada por \(f(x|\theta) = (1-\theta)^{x-1} \theta\) com \(x = 1,2,\ldots\).

A distribuição da geométrica é dada por \[p(x|\theta) = (1-\theta)^{x-1} \theta\]

O logaritmo da verosimilhança é: \[log(p(x|\mu)) = L = (x-1) log(1-\theta) + log(\theta)\]

A sua primeira derivada: \[\frac{\delta L}{\delta\theta} = \frac{1}{\theta} - \frac{(x-1)}{1-\theta}\]

E a segunda derivada: \[\frac{\delta^2 L}{\delta\theta^2} = - \frac{1}{\theta^2} - \frac{x-1}{(1-\theta)^2}\]

Assim:

\[ \begin{array}{lclr} I(\theta) & = & - E_{X|\theta} \Big[ \frac{\delta^2 L}{\delta\theta^2} \Big] & \\ & = & - E_{X|\theta} \Big[ - \frac{1}{\theta^2} - \frac{X-1}{(1-\theta)^2} \Big] & \\ & = & - \Big[ - \frac{1}{\theta^2} - \frac{\frac{1}{\theta}-1}{(1-\theta)^2} \Big] & \color{blue}{E_{X|\theta}(X) = \frac{1}{\theta}} \\ & = & \theta^{-2} (1-\theta)^{-1} & \\ \end{array} \]

A prior de Jeffreys é dada pela raiz quadrada deste valor: \[p(\theta) = \sqrt{I(\theta)} = \sqrt{\theta^{-2} (1-\theta)^{-1}} = \theta^{-1} (1-\theta)^{-1/2}\]

Exercício 24

A distribuição é dada por \[f(x|\theta) = \theta 3^\theta x^{1-\theta}\]

1. Obtenha a prior de Jeffreys:

O logaritmo da verosimilhança é: \[log(f(x|\mu)) = L = log(\theta) + log(3)\theta + (1-\theta)log(x) \]

A sua primeira derivada: \[\frac{\delta L}{\delta\theta} = \frac{1}{\theta} + log(3) - log(x)\]

E a segunda derivada: \[\frac{\delta^2 L}{\delta\theta^2} = - \frac{1}{\theta^2}\]

Assim: \[p(\theta) = \sqrt{I(\theta)} = \sqrt{- E_{X|\theta} \Big[ - \frac{1}{\theta^2} \Big]} = \frac{1}{\theta}\]

Reordenando a expressão da verosimilhança:

\[f(x|\theta) = \theta e^{log(3)x} x^{-(\theta-1)}\]

Sabemos que a verosimilhança, dado que \(x_i\) são iid:

$$

$$

Isto corresponde ao núcleo de uma distribuição gamma de parâmetros \(\alpha = n\) e \(\beta=\sum_i log(x_i) - n log(3)\).

Para uma distribuiçao gama(\(\alpha,\beta\)):

\[E(\theta|x) = \frac{\alpha}{\beta}\] \[Var(\theta|x) = \frac{\alpha}{\beta^2}\]

Exercício 25

TODO:

Exercício 26

TODO:

Exercício 27

TODO:

Exercício 28

Seja a amostra aleatória \(X_1, \ldots, X_n\) com distribuição Binomial Negativa\((r,\theta)\), \(r \gt 0, 0 \lt \theta \lt 1\), com função de probabilidade:

\[p(x|r,\theta) = \binom {x+r-1} x \theta^r (1-\theta)^x\]

com \(x = 0,1,2,\ldots\)

A conjugada da binomial negativa é a distribuiçao beta.

Assumindo uma prior uniforme, beta(a,b), temos:

\[\theta \sim beta(a,b)\]

e a regra de atualização é dada por:

\[\theta | x \sim beta(a + rn, b + \sum_i x_i)\]

com \(x = (x_1, x_2, \ldots, x_n)\).

Seja \(\theta \sim beta(2,2)\).

Para \(r=5, n=10, \sum_i x_i = 70\) temos \(\theta|x \sim beta(52,72)\).

Sejam as hipóteses \(H_0:\theta \leq 0.5\) vs. \(H_1:\theta \gt 0.5\).

Qual das hipóteses apresenta maior probabilidade à posteriori?

\[p(H_0|x) = p(\theta \leq 0.5|x) = pbeta(0.5,52,72) = 0.964\]

\[p(H_1|x) = 1 - p(H_0|x) = 0.036\]

\[O(H_1,H_0|x) = \frac{p(H_1|x)}{p(H_0|x)} = \frac{0.036}{0.964} = 0.037\]

Assumindo que as hipóteses tinham inicialmente igual probabilidade:

O factor de Bayes a favor de \(H_1\):

\[B(x) = \frac{O(H_1,H_0|x)}{O(H_1,H_0)} = \frac{0.037}{1} = 0.037\]

Não existe evidência a favor de \(H_1\).

Exercício 29

Seja \(\theta\) a probabilidade de preferir a congelação de qualidade.

Modelo de dados proposto: \(f(x|\theta) \sim binomial(\theta)\)

Não se considera as diferenças de gostos individuais dos clientes.

\(N=16,x=13\)

Para \(\theta_1 \sim beta(0.5,0.5)\), \(\theta_1 \sim beta(13.5,3.5)\)

a <- 13.5
b <- 3.5
mean <- a/(a+b)
median <- (a-1/3) / (a+b-2/3)
mode <- (a-1) / (a+b-2)
mean

## [1] 0.7941176

median

## [1] 0.8061224

mode

## [1] 0.8333333

\(p(\theta_1 \gt 0.6 | x) = 1 - p(\theta_1 \leq 0.6 | x)\)

1 - pbeta(0.6,a,b)

## [1] 0.9637777

Seja \(H_0:\theta \geq 0.6\) vs. \(H_1:\theta \lt 0.6\)

Assumimos \(p(H_0) = P(H_1) = 0.5\), logo \(O(H_0,H_1) = 1\)

\[B(x) = \frac{O(H_0,H_1|x)}{O(H_0,H_1)} = O(H_0,H_1|x) = \frac{p(\theta_1 \geq 0.6 | x)}{p(\theta_1 \lt 0.6 | x)} = \frac{0.9637}{0.036} = 26.6\]

Repetindo os cálculos para priors diferentes:

Para \(\theta_2 \sim beta(1,1)\) obtemos \(B(x) = 19\)

Para \(\theta_3 \sim beta(2,2)\) obtemos \(B(x) = 13.3\)

A prior aparenta ter um grande efeito no factor de Bayes.

Exercício 30

\[\theta \sim Gamma(1,1)\]

\[X_i \sim Poisson(\theta)\]

Para \(n=10, \sum_i x_i = 6\) temos: \[\theta | x \sim Gamma(7,11)\]

Intervalo de confiança a 90%:

qgamma(0.05,7,11)

## [1] 0.2986651

qgamma(0.95,7,11)

## [1] 1.076581

Região de credibilidade HPD a 90%:

library(TeachingDemos)

hpd(qgamma,conf=0.9, shape=7, rate=11)

## [1] 0.253163 1.005407

Exercício 31

prior: \(\theta \sim beta(1,1)\)

verosimilhança de \(x\) casos em \(n\) observações: \(x \sim binomial(\theta)\)

posterior: \(\theta | x \sim beta(x+1,n-x+1)\)

Temos as seguintes hipóteses: + \(H_0: \theta = 0.2\) + \(H_1: \theta = 0.3\)

Hipóteses à priori: + \(p(H_0) = 0.25\) + \(p(H_1) = 0.75\)

Tendo assim \(O(H_0,H_1) = \frac{0.25}{0.75} = \frac{1}{3}\)

Seja o factor de Bayes em favor de \(H_0\):

\[B(x) = \frac{O(H_0,H_1|x)}{O(H_0,H_1)}\]

Este factor pretendemos que seja maior que 1 (para favorecer \(H_0\) como é pedido no enunciado). Assim:

\[B(x) = \frac{O(H_0,H_1|x)}{O(H_0,H_1)} \gt 1 \iff \frac{O(H_0,H_1|x)}{\frac{1}{3}} \gt 1 \iff O(H_0,H_1|x) \gt \frac{1}{3}\]

Ou seja,

\[O(H_0,H_1|x) \gt \frac{1}{3} \iff \frac{p(H_0|x)}{p(H_1|x)} \gt \frac{1}{3}\iff p(H_0|x) \gt 3p(H_1|x)\]

Ora, sabendo que os dados seguem uma binomial:

• \(p(H_0|x) = p(\theta = 0.2|x) = \binom n x 0.2^x 0.8^{n-x}\)
• \(p(H_1|x) = p(\theta = 0.3|x) = \binom n x 0.3^x 0.7^{n-x}\)

Juntando tudo:

\[\binom n x 0.2^x 0.8^{n-x} \gt 3 \binom n x 0.3^x 0.7^{n-x}\]

\[0.2^x 0.8^{n-x} \gt 3 \times 0.3^x 0.7^{n-x}\]

\[x \lt \frac{n - 8.67}{4.2}\]