QuestĂŁo 01 (1o Modo)

theta <- seq(0.01,0.99,0.02)
var.est <- rep(0,length(theta))
n <- 10
licr <- ((1-2*theta)^2)*theta*(1-theta)/n
na <- 50000
for (i in 1:length(theta)) {
vec <- apply(matrix(rbinom(na*n,1,theta[i]),na,n),1,mean)
estim <- n*vec*(1-vec)/(n-1)
var.est[i] <- var(estim)
}
plot(theta,licr,type="l",ylim=c(0,0.007))
lines(theta,var.est, lty=2)

QuestĂŁo 01 (2o Modo)

VariĂąncia Exata

Seja \(Y\) uma variavel aleatória com distribuição Binomial(n,p). Utilizaremos o \(k\)-ésimo momento fatorial da distribuição binomial; \[ \begin{align} E\big[Y(Y-1)(Y-2)\cdots(Y-k+1)\big]=\frac{n!}{(n-k)!}p^k \quad (\mathbf{Prove}!) \end{align} \]

Calculemos agora alguns momentos, recursivamente:

2o Momento

\[ \begin{align} E[Y^2]=E[Y(Y-1)]=Var[Y]+(E[Y])^2=np(1-p)+(np)^2=np[(n-1)p+1] \end{align} \]

3o Momento

\[ \begin{align} \frac{n!}{(n-3)!}p^3&=E[Y(Y-1)(Y-2)]\\ &=E[Y^3]-3E[Y^2]+2E[Y]\\ \therefore \quad E[Y^3]&=\frac{n!}{(n-3)!}p^3+3E[Y^2]-2E[Y]\\ &=\frac{n!}{(n-3)!}p^3+3np[(n-1)p+1]-2np\\ &=np\big[(n-1)(n-2)p^2+3(n-1)p+1\big] \end{align} \]

4o Momento

\[ \begin{align} \frac{n!}{(n-4)!}p^4&=E[Y(Y-1)(Y-2)(Y-3)]\\ &=E[Y^4-6Y^3+11Y^2-6Y]\\ &=E[Y^4]-6E[Y^3]+11E[Y^2]-6E[Y]\\ \therefore \quad E[Y^4]&=\frac{n!}{(n-4)!}p^4+6E[Y^3]-11E[Y^2]+6E[Y]\\ &=n(n-1)(n-2)(n-3)p^4+6\big(np\big[(n-1)(n-2)p^2+3(n-1)p+1\big]\big) \\ &-11\big(np[(n-1)p+1]\big)+6np\\ &=n(n-1)(n-2)(n-3)p^4+6n(n-1)(n-2)p^3+7n(n-1)p^2+np \end{align} \]

Seja \(X_1, X_2, \cdots, X_n\) uma amostra aleatória da distribuição \(Ber(p)\). Desejamos calcular

\[ Var\bigg[\frac{n\bar{X}(1-\bar{X})}{n-1}\bigg]. \]

Utilizando o fato acima: \[ \begin{align} Var \big[T(\mathbf{X})\big] &=\bigg( \frac{n}{n-1} \bigg)^2 Var\big[\overline{X}(1-\overline{X})\big]\\ &=\bigg(\frac{n}{n-1}\bigg)^2\bigg(Var[\overline{X}]+Var[\overline{X}^2]-2Cov(\overline{X},\overline{X}^2)\bigg)\\ &=\frac{n^2}{(n-1)^2}\bigg(Var[\overline{X}]+E[\overline{X}^4]-\big(E[\overline{X}^2]\big)^2-2E\big[\big(X-E[\overline{X}]\big)\big(\overline{X}^2-E\big[\overline{X}^2\big]\bigg)\\ &=\frac{n^2}{(n-1)^2}\bigg(Var[\overline{X}]+E[\overline{X}^4]-\big(E[\overline{X}^2]\big)^2-2E\big[\overline{X}^3-\mu_2\overline{X}-\mu_1\overline{X}^2+\mu_1\mu_2\big]\bigg)\\ &=\frac{n^2}{(n-1)^2}\bigg(\frac{\mu_1(1-\mu_1)}{n}+\frac{E[Y^4]}{n^4}-\mu_2^2-\frac{2E[Y^3]}{n^3}+2\mu_2\mu_1+2\mu_1\mu_2-2\mu_1\mu_2\bigg) \end{align} \]

em que, \(\mu_1=E[\bar{X}]=p\),;; \(\mu_2=E[\bar{X}^2]\);; e;; \(Y:=\sum\limits_{i=1}^{n}X_i \sim Binomial(n,p)\).

Substituindo \(\mu_1=p\) e os resultados dos 3 momentos seguintes (Eqs. , e ) em , obtem-se, após manipulaçÔes algébricas e simplificaçÔes, a expressão seguinte:

\[ \begin{align} Var\big[T(X)\big]=(-4n^2+10n-6)p^4+(8n^2-20n+12)p^3+(-5n^2+12n-7)p^2+(n^2-2n+1)p \end{align} \] ### VariĂąncia AssintĂłtica ( Calculada pelo Aluno Pedro Brom) \[ Var\big[T(\mathbf{X}\big]=\frac{p(1-p)(4np^2-2p^2-4np+2p+n)}{(n-1)^2} \]

# Variancia exata
var.exata<-function(n,p){((-4*n^2+10*n-6)*p^4+(8*n^2-20*n +12)*p^3+
                          (-5*n^2+12*n-7)*p^2+(n^2-2*n+1)*p)/(n*(n-1)^2)}

# variancia assintotica
var.ass<-function(n,p){(((1-p)*p*(4*n*p^2-2*p^2-4*n*p+2*p+n))/(n-1)^2)}

# variancia simulada via Monte Carlo
var.sim<-function(p,N,n) {
             
            m<-length(p)
            p.res<-numeric(m)
            
         for(j in 1:m) {
               res<-numeric(N)
                for (i in 1:N){
                  X<-rbinom(n,1,p[j])
                  res[i]<-mean(X)*(1-mean(X))*(n/(n-1))
             }
        p.res[j]<-(N/(N-1))*var(res)
      }
    return(p.res)   
}

# Limite Inferior de Cramér-Rao
licr<-function(n,p){((1-2*p)^2)*p*(1-p)/n}

library(ggplot2)
ggplot()+ xlim(0,1)+
         geom_function( fun = var.sim, args = list(n=10,N=1000), aes(colour= "sim")) +
         geom_function( fun = var.ass, args = list(n=10), aes(colour= "ass")) +
         geom_function( fun = var.exata, args = list(n=10), aes(colour= "exata"))+
         geom_function( fun = licr, args = list(n=10), aes(colour= "licr"))+
         labs(title="VariĂąncia", x=bquote(p ~ (n==10)), y="Valores")+
  theme(plot.title = element_text(size=12,hjust=0.5))

QuestĂŁo 03

n<-1000 ; m<-10 ; p<-0.4
x<-rbinom(n,m,p)

licr<-function(p,n,m){ p*(1-p)/(n*m)}
varian<-licr

library(ggplot2)
ggplot()+
  xlim(0.1,0.90)+
  geom_function(fun = licr, args = list(n=1000,m=10))+
  labs(title="Limite Inferior de Cramér-Rao",  x="p",  y="LICR")+
  theme(plot.title = element_text(h=0.5))

QuestĂŁo 04

licr<-function(n,lambda){log(2)^2/(n*lambda^2)}

ggplot()+
  xlim(0.1,1)+
  geom_function( fun = licr, args = list(n=10),colour= "red")+
  labs(title="Lmite Inferior de Cramér-Rao", 
       x=bquote(lambda ~ (n==10)), y="LICR")+
  theme(plot.title = element_text(size=12,hjust=0.5))

QuestĂŁo 05

## 1o Modo ( R base)
theta <- seq(-3,3,0.01)
n <- 10
var.est <- pnorm(theta,0,1)*(1-pnorm(theta,0,1))/n
licr <- exp(-theta^2)/(2*pi*n)
plot(theta,licr,type="l",ylim=c(0,0.03))
lines(theta,var.est,lty=2)

## 2o Modo ( Com ggplot do R )
library(ggplot2)
ggplot() + xlim(-3,3) +
  geom_function(fun = ~ exp(- (.x)^2 )/(2*pi*n), args=list(n=10),aes(colour = "LICR") ) +
  geom_function(fun = ~ pnorm(.x)*(1-pnorm(.x))/n, args=list(n=10),aes(colour = "Var") ) +
  labs(title="VariĂąncia X LICR", x=bquote(theta), y="Valores") +
  theme(plot.title = element_text(h=0.5))

QuestĂŁo 06

Para uma amostra aleatĂłria de tamanho \(n\) da dist. Poi(\(\lambda\)).

  1. Encontre o estimador nĂŁo tendencioso de mĂ­nima variĂąncia de \(\tau(\lambda)=e^{-\lambda}\).

  2. Obtenha a variùncia deste estimador e compare com o LICR. Compare (a variùncia obtida) também com a variùncia do EMV (que é tendencioso).

Solução

Item a) JĂĄ feito pelo professor.

Item b) Observe que \(LICR=\frac{\lambda e^{-2\lambda}}{n}\) (calcule!). E, considerando uma amostra aleatória \(X_1,\ldots,X_n\) da distribuição \(Pois(\lambda)\), temos que \[S:=\sum_{i=1}^{n}X_i \sim Pois(n\lambda)\]

Aplicaremos o seguinte fato:

Para todo \(a\neq 0\), vale \[ \begin{align} E(a^S)&=\sum_{s=0}^{\infty}a^sP(S=s)\\ &=\sum_{s=0}^{\infty}a^s \frac{e^{-n\lambda}(n\lambda)^s}{s!}\\ &=e^{-n\lambda}\sum_{s=0}^{\infty}\frac{(n\lambda a)^s}{s!}\\ &=e^{-n\lambda}e^{n\lambda a}=e^{n\lambda(a-1)} \end{align} \]

Variùncia do Estimador de Måximo Verossimilhança (\(EMV=e^{-\bar{X}}\))

Para obter a variĂąncia, tomemos \(a=e^{-1/n}\) para encontrar \(E(e^{-\bar{X}})\) e \(a=e^{-2/n}\) para calcular \(E(e^{-2\bar{X}})\).

Com isso,

\[ \begin{align} Var\big[e^{-\bar{X}}\big]&=E\big[e^{-2\bar{X}}\big] - \big\{E\big[e^{-\bar{X}}\big]\big\}^2\\ &=E\big[e^{-\frac{2}{n}S}\big] -\big\{E\big[e^{-\frac{1}{n}S}\big]\big\}^2\\ &=e^{n\lambda\big(e^{-2/n}-1\big)} - e^{2n\lambda\big(e^{-1/n}-1\big)} \qquad \mbox{Prof.,lĂĄ na lista, hĂĄ erros nos expoentes!} \end{align} \]

VariĂąncia do Estimador \(T(\mathbf{X})=\big(\frac{n-1}{n}\big)^S\)

Utilizemos \(a=\frac{n-1}{n}\) para encontrar a esperança e \(a=\big(\frac{n-1}{n}\big)^2\) para o 2o momento, conforme abaixo

\[ \begin{align} E\bigg[\bigg(\frac{n-1}{n}\bigg)^{S}\bigg]&=e^{n\lambda\big(\frac{n-1}{n}-1\big)}=e^{-\lambda}\\\\ E\bigg[\bigg(\frac{n-1}{n}\bigg)^{2S}\bigg]&=e^{n\lambda\big[ \big(\frac{n-1}{n}\big)^2-1\big]}\\ &=e^{-2\lambda \;+\;\lambda/n}\\\\ \therefore \quad Var\bigg[\bigg(\frac{n-1}{n}\bigg)^{S}\bigg]&=e^{-2\lambda\;+\;\lambda/n}\;-\;\big(e^{-\lambda}\big)^2\\ &=e^{-2\lambda}(e^{\lambda/n}-1) \end{align} \]

Comparando VariĂąncia do Estimador T com o LICR.

\[ R=\frac{Var(T)}{LICR}=\frac{e^{-2\lambda}(e^{\lambda/n}-1)}{\lambda e^{-2\lambda}/n}=\frac{n(e^{\lambda/n}-1)}{\lambda} \]

Logo,

\[ R\geq 1, \quad (i.e, Var(T)\geq LICR) \iff e^{\lambda/n}-1\geq \lambda/n \]

Mostraremos que a desigualdade \(e^{\lambda/n}-1\geq \lambda/n\) Ă© vĂĄlida para todo \(\frac{\lambda}{n}>0\).

1a Demonstração: Seja \(x=\lambda/n\). Como \(\lambda>0\) e \(n>0\), então \(x>0\). Uma vez que \(f(t)=e^{t}\) é uma função contínua crescente, temos que \(e^{x}\)>1. Logo, pelo Teorema do Valor Médio, existe \(h\in (0,x)\) tal que, para todo \(x>0\), \[ \begin{align} e^{x}-e^0&=e^{h}(x-0)\\ \therefore \quad e^{x}-1&=e^{h}x \iff e^{x}=1+e^hx > 1+x \quad \because \quad e^h>0.\\ \end{align} \] 2a Demonstração: Seja \(x=\lambda/n>0\). Como \(f(t)=e^{t}\) é uma função contínua crescente, segue:

\[ \begin{align} \forall\; x>0, \quad 0=\int_{0}^{x}0dt &\leq \int_{0}^{x}(e^t-1)dt\\ &=e^{t}\bigg|_{0}^{x}- t\big|_{0}^{x}\\ &=e^{x}-x-1.\\ \therefore \quad & e^{-x}-1\geq x. \end{align} \]

## Limite Inferior de Cramér-Rao
licr<-function(n,lambda){ lambda*exp(-2*lambda)/n}

## VariĂąncia de estimador T (est de tau(lambda))

var.T<-function(n,lambda){ exp(-2*lambda)*(exp(lambda/n) - 1)}

## Comparando a variancia do estimador T com o LICR.
 # Importante: nao comparamos LICR com a var do emv, por este ser um est. viesado.

library(ggplot2)
ggplot() + xlim(0.1,5) +
  geom_function( fun = var.T, args = list(n=10), aes(colour= "Est.T")) +
  geom_function( fun = licr, args = list(n=10), aes(colour= "LICR")) +
  labs(title="VariĂąncia X LiCR", x=bquote(lambda~(n==10)), y="Valores") +
  theme(plot.title = element_text(h=0.5))

Comparando as variĂąncias dos estimadores T e MV.

## VariĂąncia dos estimadores

var.T<-function(n,lambda){ exp(-2*lambda)*(exp(lambda/n) - 1)}

var.mv<-function(n,lambda){ exp(n*lambda*(exp(-2/n)-1)) - exp(2*n*lambda*(exp(-1/n)-1))}

library(ggplot2)
ggplot() + xlim(0.1,5) +
  geom_function( fun = var.T, args = list(n=10), aes(colour= "Est.T")) +
  geom_function( fun = var.mv, args = list(n=10), aes(colour= "Est.MV")) +
  labs(title="VariĂąncia", x=bquote(lambda~(n==10)), y="Var") +
  theme(plot.title = element_text(h=0.5))

Comparando os EQMs dos estimadores T e MV.

Esp.mv<-function(n,lambda){exp(n*lambda*(exp(-1/n)-1))}
EQM.mv<-function(n,lambda){var.mv(n,lambda) + ( Esp.mv(n,lambda) - exp(-lambda) )^2}

Esp.T<-function(n,lambda){exp(-lambda)}
EQM.T<-function(n,lambda){var.T(n,lambda) + ( Esp.T(n,lambda) - exp(-lambda) )^2}

## Comparando  os EQMs

library(ggplot2)
ggplot() + xlim(0.1,5) +
  geom_function( fun = EQM.T, args = list(n=10), aes(colour= "Est.T")) +
  geom_function( fun = EQM.mv, args = list(n=10), aes(colour= "Est.MV")) +
  labs(title="EQM", x=bquote(lambda~(n==10)), y="eqm") +
  theme(plot.title = element_text(h=0.5))

## 1o Modo
theta <- 2
n <- 10
x <- rbeta(n,theta,1)
theta.vec <- seq(0.01,5,0.01)
veros <- rep(0,length(theta.vec))
for (i in 1:length(theta.vec))
veros[i] <- prod(dbeta(x,theta.vec[i],1))
plot(theta.vec,veros,type="l")
lambda <- 2
r.par <- 2
priori <- dgamma(theta.vec,r.par,rate=lambda)
lines(theta.vec,priori,lty=2)
posteriori <- dgamma(theta.vec,(n+r.par),rate=(lambda-sum(log(x))))
lines(theta.vec,posteriori,lty=3)
legend(3,2,legend=c("Veross.","priori","posteriori"),lty=c(1,2,3))

## 2o Modo