Enunciado dos Exercícios

Exercício 1

De acordo com o enunciado temos os seguintes dados:

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:

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)\)

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:

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

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)

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