Estatística Computacional II - Lista 4

Aluno: Rafael Cabral Fernandez

Introdução

Este documento refere-se a lista 4 de exercícios proposta pelo professor Gustavo Rocha, da cadeira de Estatística Computacional 2, como cômputo de ensino da disciplina.

1

Para efetuar o algoritmo EM, precisamos primeiro avaliar o valor esperado das variáveis latentes em função das variáveis conhecidas, isto é:

\(Q(\theta,\theta^*) = E_{z|x,\theta}\bigg[ln(f(X,Z|\theta))\bigg]=E_{z|x,\theta}\bigg[ln\bigg(\cfrac{197}{z_1!z_2!18!20!34!}\bigg(\cfrac{1}{2}\bigg)^{z_1}\bigg(\cfrac{\theta}{4}\bigg)^{z_2}\bigg(\cfrac{1}{4}-\cfrac{\theta}{4}\bigg)^{18}\bigg(\cfrac{1}{4}-\cfrac{\theta}{4}\bigg)^{20}\bigg(\cfrac{1}{2}\bigg)^{34} \bigg) \bigg]\)

\(\propto\)

\(E_{z|x,\theta}\bigg[z_1ln(\cfrac{1}{2}) + z_2ln(\cfrac{\theta}{4}) + 18ln(\cfrac{1}{4}-\cfrac{\theta}{4}) +20ln(\cfrac{1}{4}-\cfrac{\theta}{4}) + 34ln(\cfrac{\theta}{4})\bigg]\)

=

\(E_{z|x,\theta}[z_1]ln(\cfrac{1}{2}) + E_{z|x,\theta}[z_2]ln(\cfrac{\theta}{4}) + 18ln(\cfrac{1}{4}-\cfrac{\theta}{4}) +20ln(\cfrac{1}{4}-\cfrac{\theta}{4}) + 34ln(\cfrac{\theta}{4})\)

Faz-se necessário definir as esperanças descritas acima. Sabe-se que \(Z_2|X_2,X_3,X_4 \sim bin(125,\cfrac{\theta^*}{2+\theta^*})\)

Logo, para efetuar a segunda etapa, devemos maximizar.

\(\cfrac{\partial Q(\theta,\theta^*)}{\partial \theta} = \cfrac{a}{\theta} - \cfrac{18}{1-\theta} - \cfrac{20}{1-\theta} + \cfrac{34}{\theta} = 0\)

\(=>\)

\(\hat{\theta} = \cfrac{34+a}{a+72}\), onde \(a = 125\bigg(\cfrac{\theta^*}{2+\theta^*}\bigg)\)

theta <- vector()

theta[1] <- 0.5
theta[2] <- (125*theta[1]/(2+theta[1]) + 34)/(125*theta[1]/(2+theta[1])+72)

i <- 2

while (abs(theta[i]-theta[(i-1)])>10^(-20)) {
  theta[i+1] <- (125*theta[1]/(2+theta[1]) + 34)/(125*theta[1]/(2+theta[1])+72)
  i <- i + 1
}

theta

## [1] 0.5000000 0.6082474 0.6082474

2

dados_lista4_questao2 <- read.table("C:/Users/Rafael/Desktop/dados_lista4_questao2.txt", quote="\"", comment.char="")

x = dados_lista4_questao2$V1

beta1 = function(x){
  196/mean(x)
}

beta2 = function(x){
  mean(x)/var(x)
}

b1 = mean(length(x)*196/mean(x)-(length(x)-1)*jackknife(x,beta1)$jack.values)
b2 = mean(length(x)*mean(x)/var(x)-(length(x)-1)*jackknife(x,beta2)$jack.values)

varb1 = jackknife(x,beta1)$jack.se
varb2 = jackknife(x,beta2)$jack.se

b1;b2;varb1;varb2

## [1] 28.03468

## [1] 28.36977

## [1] 0.06289077

## [1] 1.230796

O estimador \(b_1\) possui vício e variância menor, por conseguinte, também deterá um menor Erro Quadrático Médio, sendo, portanto, um estimador mais preciso.

b1 = mean(bootstrap(x,10000,beta1)$thetastar)
b2 = mean(bootstrap(x,10000,beta2)$thetastar)

intervalob1 = quantile(bootstrap(x,10000,beta1)$thetastar,c(.025,0.975))
intervalob2 = quantile(bootstrap(x,10000,beta2)$thetastar,c(.025,0.975))

varb1 = var(bootstrap(x,10000,beta1)$thetastar)
varb2 = var(bootstrap(x,10000,beta2)$thetastar)


b1;varb1;intervalob1;

## [1] 28.03367

## [1] 0.003942984

##    2.5%   97.5% 
## 27.9099 28.1594

b2;varb2;intervalob2;

## [1] 28.52873

## [1] 1.540794

##     2.5%    97.5% 
## 26.19652 31.06150

beta = 196/mean(x)

amostra = rgamma(10000,196,beta)

b1 = mean(bootstrap(amostra,10000,beta1)$thetastar)
b2 = mean(bootstrap(amostra,10000,beta2)$thetastar)

intervalob1 = quantile(bootstrap(amostra,10000,beta1)$thetastar,c(.025,0.975))
intervalob2 = quantile(bootstrap(amostra,10000,beta2)$thetastar,c(.025,0.975))


b1;intervalob1;

## [1] 28.06338

##     2.5%    97.5% 
## 28.02467 28.10442

b2;intervalob2

## [1] 27.77029

##     2.5%    97.5% 
## 27.00520 28.54997

hist(amostra,breaks=50,probability = T)
curve(dgamma(x,196,b1),0,10,add=T,col="red")
curve(dgamma(x,196,b2),0,10,add=T,col="blue")

A curva em vermelho é a obtida com a estimativa de b1, a em azul, com a estimativa de b2. Como é possível observar visualmente, a curva da densidade Gama com o parâmetro obtido atráves do estimador b1 é a que mais se aproxima dos dados reais.

3

Parar estimar o \(p\), faremos \(p^j = \sum_{i=1}^{1000}\cfrac{w(x_i)p^{j-1}}{w(x_i)p^{j-1} + g(x_i)(1-p)^{j-1}}\)

Onde,

\(w(x)\) é a densidade de uma distribuição Gama com parâmetros (4,\(\beta_1\)) e \(g(x)\) é a densidade de uma distribuição Gama com parâmetros (196,\(\beta_2\))

No entanto, será necessário estimar \(\theta=(p,\beta_1,\beta_2)\).

Para estimação de \(\beta_1\), faça:

\(\beta_1^{(j)} = \cfrac{4\sum_{i=1}^{1000}a_i^{(j)}}{\sum_{i=1}^{1000}a_i^{(j)}x_i}\)

E para estimação de \(\beta_2\), faça:

\(\beta_2^{(j)} = \cfrac{196\sum_{i=1}^{1000}(1-a_i^{(j)})}{\sum_{i=1}^{1000}(1-a_i^{(j)}) x_i}\)

x <- read.table("C:/Users/Rafael/Desktop/dados_lista4_questao3.txt", quote="\"", comment.char="")
x <- x[,1]

a=c()
theta0=c(0.5,1,30)
theta=c()
j=0

norma_r=Inf
while(norma_r>10^(-20)){
  for(i in 1:length(x)){
    pw=theta0[1]*dgamma(x[i],4,theta0[2])
    a[i]=pw/(pw+(1-theta0[1])*dgamma(x[i],196,theta0[3]))}
  theta[1]=sum(a)/length(x)
  theta[2]=(4*sum(a))/(sum(a*x))
  theta[3]=196*(sum(1-a)/sum(x*(1-a)))
  j=j+1
  r=theta-theta0
  norma_r=sqrt(sum(r^2))
  theta0=theta}

  
print(theta)

## [1]  0.7251579  1.9975157 28.1315973

print(j)

## [1] 20

hist(x,breaks=50,probability = T)
curve(theta[1]*dgamma(x,4,theta[2])+(1-theta[1])*dgamma(x,196,theta[3]),0,10,add=T)

Como a curva obtida através através da densidade com os parâmetros estimados se aproxima do histograma da amostra, pode-se concluir que a qualidade da estimativa é boa.