Lista 2 Inferência - Computacional
Adolfo M. D. Silva
21/03/2021
- Questão 01 (1o Modo)
- Questão 01 (2o Modo)
- Questão 03
- Questão 04
- Questão 05
- Questão 06
- Variância do Estimador de Máximo Verossimilhança (\(EMV=e^{-\bar{X}}\))
- Variância do Estimador \(T(\mathbf{X})=\big(\frac{n-1}{n}\big)^S\)
- Comparando Variância do Estimador T com o LICR.
- Comparando as variâncias dos estimadores T e MV.
- Comparando os EQMs dos estimadores T e MV.
- Estimador UMVU
- Outro estimador não viesado de \(\tau(\lambda)=(1+\lambda)e^{-\lambda}\)
- Questão 07d)
- Questão 08
- Questão 09
- Questão 13
- Questão 16
- Questão 18
- Questão 19
- Questão 20
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 os primeiros 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]. \]
Solução: Defina-se \(T(\mathbf{X})= \frac{n\bar{X}(1-\bar{X})}{n-1}\). Logo, \[ \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\bigg), \end{align} \]
em que \(\mu_1=E[\bar{X}]=p \quad {,}\quad \mu_2=E[\bar{X}^2]=\frac{p(1-p)}{n}+p^2 \quad \mbox{e} \quad Y:=\sum\limits_{i=1}^{n}X_i \sim Binomial(n,p)\).
Substituindo \(\mu_1=p\) e os resultados, referentes aos 2o, 3o e 4o momentos, na útlima equação, obtem-se, após manipulações algébricas e simplificações, a fórmula seguinte: \[ \begin{align} Var\big[T(\mathbf{X})\big]=\frac{(-4n^2+10n-6)p^4+(8n^2-20n+12)p^3+(-5n^2+12n-7)p^2+(n^2-2n+1)p}{n(n-1)^2}. \end{align} \]
Variância Assintótica (Aproximação com uso da dist. Normal)
( Esta variância foi 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}
# Graficos
library(ggplot2)
g1<-ggplot()+ xlim(0,1)+
geom_function( fun = var.sim, args = list(n=10,N=1000), aes(colour= "simulada")) +
geom_function( fun = var.ass, args = list(n=10), aes(colour= "assintótica")) +
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))
g2<-ggplot()+ xlim(0,1)+
geom_function( fun = var.sim, args = list(n=30,N=1000), aes(colour= "simulada")) +
geom_function( fun = var.ass, args = list(n=30), aes(colour= "assintótica")) +
geom_function( fun = var.exata, args = list(n=30), aes(colour= "exata"))+
geom_function( fun = licr, args = list(n=30), aes(colour= "LICR"))+
labs(title="Variância", x=bquote(p ~ (n==30)), y="Valores")+
theme(plot.title = element_text(size=12,hjust=0.5))
library(ggpubr)
ggarrange(g1, g2, ncol=2, nrow=1, common.legend = TRUE, legend="bottom")
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="Limite 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
Seja \(X_1,X_2,\cdots,X_n\) uma amostra aleatória da distribuição \(N(\theta,1)\).
Desejamos calcular o LICR para a variância do estimador não viesado de \(P(X>0)\).
Solução: Defina-se \(Z:=X-\theta\) e \(P(Z\leq \theta)=\Phi(\theta)\). Desse modo, para todo \(\theta\) real, temos \[ \begin{align} P(X>0)&=P\bigg(\frac{X-\theta}{\sqrt{1}}>\frac{0-\theta}{\sqrt{1}}\bigg)\\ &=P(Z>-\theta)\\ &=P(Z\leq \theta) \quad \because \quad (\mbox{simetria da dist. Normal})\\ &=\Phi(\theta). \end{align} \] Portanto, \[ \begin{align} LICR=\frac{\big[\Phi^{'}(\theta)\big]^2}{nI_{F}(\theta)}=\frac{\big[\phi_N(\theta)\big]^2}{n\times1}=\frac{\bigg(\frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}\theta^2}\bigg)^2}{n}=\frac{1}{2\pi n }e^{-\theta^2}, \quad \theta \in (-\infty,\infty). \end{align} \]
Obs.: ( Façam os outros itens! )
## 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\)).
Encontre o estimador não tendencioso de mínima variância de \(\tau(\lambda)=e^{-\lambda}\).
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) Sabemos 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)\), tem-se \[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: De fato. 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>1.\\ \end{align} \] 2a Demonstração: Com efeito, seja \(x=\lambda/n>0\). Como \(f(t)=e^{t}\) é uma função contínua crescente, segue que
\[ \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))
c) Encontrando o estimador UMVU
Solução
Estimador UMVU
Seja \(S:=\sum\limits_{j=1}^{n}X_j\) e \(s\in\{0,1,\ldots\}\) os valores assumidos por \(S\).
Logo, \(S:=\sum\limits_{j=1}^{n}X_j \sim Poi(n\lambda)\). Desse modo,
\[ \begin{align} E\bigg[I\big\{X_i=0\big\}+I\big\{X_i=1\big\}\bigg|S=s\bigg]&=E\bigg[I\big\{X_i=0\big\}\bigg|\sum\limits_{j=1}^{n}X_i=s\bigg]+E\bigg[I\big\{X_i=1\big\}\bigg|\sum\limits_{j=1}^{n}X_i=s\bigg]\\ & = \frac{P\big(X_i=0; \sum\limits_{j=1}^{n}X_j=s\big)}{P\big(\sum\limits_{j=1}^{n}X_j=s\big)} + \frac{P\big(X_i=1; \sum\limits_{j=1}^{n}X_j=s\big)}{P\big(\sum\limits_{j=1}^{n}X_j=s\big)}\\ &=\frac{P\big(X_i=0\big)P\big(\sum\limits_{j\neq i}X_j=s\big)}{P\big(\sum\limits_{j=1}^{n}X_j=s\big)}+ \frac{P\big(X_i=1\big)P\big(\sum\limits_{j\neq i}X_j=s-1\big)}{P\big(\sum\limits_{j=1}^{n}X_j=s\big)}\\ &=\frac{e^{-\lambda}e^{-\lambda(n-1)}[\lambda(n-1)]^s/s!}{e^{-n\lambda}[n\lambda]^s/s!}+ \frac{\lambda e^{-\lambda}e^{-\lambda(n-1)}[\lambda(n-1)]^{s-1}/(s-1)!}{e^{-n\lambda}[n\lambda]^s/s!} \quad \because \quad \sum\limits_{j\neq i}X_j\sim Poi\big(\lambda(n-1)\big)\\ &=\bigg(\frac{n-1}{n}\bigg)^s + \frac{s}{n}\bigg(\frac{n-1}{n}\bigg)^{s-1}. \end{align} \]
Portanto, \[ R(\mathbf{X})=\widehat{\tau(\lambda)}_{umvu}=\bigg(\frac{n-1}{n}\bigg)^S + \frac{S}{n}\bigg(\frac{n-1}{n}\bigg)^{S-1}. \]
Variância do Estimador UMVU
Calculemos apenas o 2o momento. ( Verifique (confirme) que R é estimador não viesado de \((1+\lambda)e^{-\lambda}\))
\[ \begin{align} E\big\{\big[R(\mathbf{X})\big]^2\big\}&=E\bigg[\bigg(\frac{n-1}{n}\bigg)^S + \frac{S}{n}\bigg(\frac{n-1}{n}\bigg)^{S-1}\bigg]^2\bigg\}\\ &=E\bigg[\bigg(\frac{n-1}{n}\bigg)^{2S}+\frac{2S}{n}\bigg(\frac{n-1}{n}\bigg)^{2S-1}+\frac{S^2}{n^2}\bigg(\frac{n-1}{n}\bigg)^{2S-2}\bigg]\\ &=E\bigg[\bigg(\frac{n-1}{n}\bigg)^{2S}\bigg]+\frac{2}{n}\bigg(\frac{n-1}{n}\bigg)^{-1}E\bigg[\bigg(\frac{n-1}{n}\bigg)^{2S}S\bigg]+ \frac{1}{n^2}\bigg(\frac{n-1}{n}\bigg)^{-2}E\bigg[\bigg(\frac{n-1}{n}\bigg)^{2S}S^2\bigg] \end{align}\\ \]
Mas, \[ \begin{align} E\bigg[\bigg(\frac{n-1}{n}\bigg)^{2S}S\bigg]&=\sum\limits_{s=0}^{\infty}\bigg(\frac{n-1}{n}\bigg)^{2s}s \;e^{-n\lambda}\frac{(n\lambda)^s}{s!}\\ &=(n\lambda)e^{-n\lambda}\bigg(\frac{n-1}{n}\bigg)^2\sum\limits_{s=1}^{\infty}\bigg\{\bigg[\bigg(\frac{n-1}{n}\bigg)^2\bigg]^{s-1}\frac{(n\lambda)^{s-1}}{(s-1)!}\bigg\}\\ &=\lambda e^{-n\lambda}\frac{(n-1)^2}{n}\sum\limits_{t=0}^{\infty}\frac{\big[\frac{\lambda(n-1)^2}{n}\big]^t}{t!}\\ &=\lambda e^{-n\lambda}\frac{(n-1)^2}{n}e^{\frac{\lambda(n-1)^2}{n}}\\ &=\frac{(n-1)^2}{n}\lambda e^{-2\lambda+\lambda/n}, \end{align} \]
\[ \begin{align} E\bigg[\bigg(\frac{n-1}{n}\bigg)^{2S}S^2\bigg]= \frac{\lambda\;(n-1)^2\;e^{\lambda(1/n-2)}\big[\lambda(n-1)^2+n\big]}{n^2} \quad (\mathbf{Prove!}), \end{align} \] e, por fim \[ E\bigg[\bigg(\frac{n-1}{n}\bigg)^{2S}\bigg] = e^{n\lambda\bigg[\bigg(\frac{n-1}{n}\bigg)^2-1\bigg]}=e^{-2\lambda+\lambda/n}. \] Logo, \[ \begin{align} E\big\{\big[R(\mathbf{X})\big]^2\big\}&=E\bigg[\bigg(\frac{n-1}{n}\bigg)^{2S}\bigg]+\frac{2}{n}\bigg(\frac{n-1}{n}\bigg)^{-1}E\bigg[\bigg(\frac{n-1}{n}\bigg)^{2S}S\bigg]+\frac{1}{n^2}\bigg(\frac{n-1}{n}\bigg)^{-2}E\bigg[\bigg(\frac{n-1}{n}\bigg)^{2S}S^2\bigg]\\ &=e^{-2\lambda+\lambda/n}+\frac{2}{n}\bigg(\frac{n}{n-1}\bigg)\frac{(n-1)^2}{n}\lambda e^{-2\lambda+\lambda/n}+\frac{1}{n^2}\bigg(\frac{n}{n-1}\bigg)^2\frac{\lambda\;(n-1)^2\;e^{\lambda(1/n-2)}\big[\lambda(n-1)^2+n\big]}{n^2}\\ &=\frac{e^{\lambda(1/n-2)}\big[\lambda^2(n-1)^2+n^2+\lambda n(2n-1)\big]}{n^2}. \end{align} \]
Portanto, \[ \begin{align} Var\big[R(\mathbf{X})\big]&=\frac{e^{\lambda(1/n-2)}\big[\lambda^2(n-1)^2+n^2+\lambda n(2n-1)\big]}{n^2} - \big[(1+\lambda)e^{-\lambda}\big]^2 \\ &=e^{-2\lambda}\big\{e^{\lambda/n}\big[\lambda^2(1-1/n)^2+1+\lambda(2-1/n)\big] - (1+\lambda)^2\big\}. \end{align} \]
LICR para variância dos estimadores não viesados de \(\tau(\lambda)=(1+\lambda)e^{-\lambda}\)
\[ LICR=\lambda^3e^{-2\lambda}/n. \]
Outro estimador não viesado de \(\tau(\lambda)=(1+\lambda)e^{-\lambda}\)
\[ \widehat{\tau(\lambda)}=\frac{1}{n}\sum\limits_{i=1}^{n}\bigg(I\big\{X_i=0\big\}+I\big\{X_i=1\big\}\bigg). \]
Variância deste estimador: (já calculada)
\[ \begin{align} Var\big[\widehat{\tau(\lambda)}\big]=\frac{(1+\lambda)e^{-\lambda}\big[1-(1+\lambda)e^{-\lambda}\big]}{n}. \end{align} \]
Implementação Computacional Item e)
## Variancia do outro est. nao viesado de (1+lambda)e^(-lambda).
var.est.tau<-function(n,l){(1+l)*exp(-l)*(1-(1+l)*exp(-l))/n}
## Variancia do estimador UMVU (estimador R(X))
var.est.R<-function(n,l){ exp(-2*l) * ( exp(l/n)* ( 1+l*(2-1/n) + l^2*(1-1/n)^2 ) - (1+l)^2 ) }
## Simulacao Monte Carlo da variancia do UMVU
var.sim<-function(N,n,l){
p.res<-numeric(length(l))
for (j in 1:length(l)) {
res<-numeric(N)
for (i in 1:N){
X<-rpois(n,l[j])
res[i]<-((n-1)/n)^(n*mean(X)) + mean(X)*( ((n-1)/n)^(n*mean(X)-1) )
}
p.res[j]<-var(res)
}
return(p.res)
}
## Limite Inferior de cramér-Rao
licr<-function(n,l){l^3*exp(-2*l)/n}
## Graficos
library(ggplot2)
g1<-ggplot() + xlim(0,7)+
geom_function( fun = var.est.tau, args = list(n=5), aes(colour= "Var. est.tau")) +
geom_function( fun = var.sim, args = list(n=5,N=1000), aes(colour="Var. Simulada")) +
geom_function( fun = licr, args = list(n=5), aes(colour= "LICR"))+
geom_function( fun = var.est.R, args = list(n=5), aes(colour = "Var. UMVU"))+
labs(title="Variância X LICR", x=bquote(lambda ~ (n==5~e~N==1000)), y="Valores")+
theme(plot.title = element_text(size=12,hjust=0.5))
g2<-ggplot() + xlim(0,7)+
geom_function( fun = var.est.tau, args = list(n=5), aes(colour= "Var. est.tau")) +
geom_function( fun = var.sim, args = list(n=5,N=1e4), aes(colour= "Var. Simulada")) +
geom_function( fun = licr, args = list(n=5), aes(colour= "LICR"))+
geom_function( fun = var.est.R, args = list(n=5), aes(colour = "Var. UMVU"))+
labs(title="Variância X LICR", x=bquote(p ~ (n==5~e~N==10000)), y="Valores")+
theme(plot.title = element_text(size=12,hjust=0.5))
library(ggpubr)
ggarrange(g1, g2, ncol=2, nrow=1, common.legend = TRUE, legend="bottom")
Questão 07d)
Por definição, \[ f_{X_1}(x)=\begin{cases} \frac{2x}{\theta^2}, & x\in (0,\theta)\\ 0, & c.c \end{cases} \]
Integrando esta função de densidade de probabilidade, obtem-se \[ F_{X_1}(x)=\begin{cases} 0, & x<0,\\ \frac{x^2}{\theta^2}, & 0\leq x<\theta,\\ 1, & x\geq \theta. \end{cases} \]
Sendo \(F_{X_{(n)}}(x)\overset{(i.i.d)}{=}\big[F_{X_{1}}(x)\big]^n\), tem-se \(f_{X_{(n)}}(x)=\frac{d}{dx}F_{X_{(n)}}(x)=n\big[F_{X_{1}}(x)\big]^{n-1}f_{X_1}(x)\).
Daí, \[ f_{X_{(n)}}(x)=\begin{cases} \frac{2n}{\theta^{2n}}x^{2n-1}, & x\in (0,\theta),\\ 0, & c.c. \end{cases} \]
Defina-se \(Y:=-log\bigg(\dfrac{X_{(n)}}{\theta}\bigg)\).
Seja \(x\in(0,\theta)\) o valor assumido pela variável aleatória \(X_{(n)}\). Desse modo,
\(y=-log\big(\frac{x}{\theta}\big)\) com \(y\in (0,\infty)\) e \(x=\theta e^{-y}\). Além disso,
\[ \frac{d}{dy}x=\frac{d}{dy}(\theta e^{-y})=-\theta e^{-y}, \;\;\; y>0 \; \wedge \; \theta >0. \]
Pela Transformação Jacobiana, temos: \[ \begin{align} f_Y(y)&=f_{X_{(n)}}(x)\bigg|\frac{d}{dy}x \bigg|\\ &=f_{X_{(n)}}(\theta e^{-y})\;\big|-\theta e^{-y}\big|\\ &=\frac{2n}{\theta^{2n}}(\theta e^{-y})^{2n-1}\theta e^{-y}\\ &=2ne^{-2ny}, \quad y>0. \end{align} \]
Logo, \[ E\bigg[-log\bigg(\frac{X_{(n)}}{\theta}\bigg)\bigg]=E\big[Y\big]=\int\limits_{0}^{\infty}y\;\underbrace{2ne^{-2ny}}_{exp(2n)}\;dy=\frac{1}{2n}. \]
Portanto, \[ E\big[-\big(\;log(X_{(n)})-log(\theta)\;\big)\big]=\frac{1}{2n} \iff E\big[log(X_{(n)})]=log(\theta)-\frac{1}{2n}. \]
E, por fim, uma vez que o máximo da amostra é uma estatística suficiente e completa (mostre!), segue que o estimador UMVU para \(log(\theta)\) é \[ \widehat{log(\theta)}_{umvu}=log\big(X_{(n)}\big)+\frac{1}{2n}. \]
# Observacao: Geracao de amostra (pela inversa da FDA)
set.seed(12345)
v.size<-c(10,100,500,1e3,1e4,1e5) # vetor de tamanhos das amostras
N<-1000 # repeticoes
theta<-3 # parametro
est.umvu<-numeric(length(v.size)) # vetor com as ests umvu
for (j in 1:length(v.size)){
x_max<-numeric(N)
for (i in 1:N){
u<-runif(v.size[j],0,1)
x<-theta*sqrt(u) # amostra gerada (observada)
x_max[i]<-max(x)
}
est.umvu[j]<-mean(log(x_max)+1/(2*v.size[j])) # estimativa umvu
}
## valor do parametro
(parametro<-log(3))## [1] 1.098612## estimativas
(df<-data.frame(v.size,est.umvu))## v.size est.umvu
## 1 1e+01 1.097288
## 2 1e+02 1.098813
## 3 5e+02 1.098624
## 4 1e+03 1.098629
## 5 1e+04 1.098611
## 6 1e+05 1.098612Variância de \(log\big(X_{(n)}\big)\)
\[ \begin{align} M_{log(X_{(n)})}(t)&:=E\bigg[e^{t\;log(X_{(n)})}\bigg]\\ &=E\big[X_{(n)}^t\big]\\ &=\int\limits_{0}^{\theta}x^tf_{X_{(n)}}dx\\ &=\int\limits_{0}^{\theta}x^t\frac{2n}{\theta^{2n}}x^{2n-1}dx\\ &=\frac{2n}{2n+t}\theta^t, \qquad \theta>0. \end{align} \]
calculando as 1a e 2a derivadas da função momento em relação a \(t\) e avaliando-as no tempo \(t=0\), obtem-se
\[ E\big[log(X_{(n)}) \big]= \frac{d}{dt}M_{log(X_{(n)})}(t)\bigg|_{t=0}=-\frac{1}{2n}+log(\theta) \]
e
\[ E\big[log(X_{(n)})^2 \big]= \frac{d^2}{dt^2}M_{log(X_{(n)})}(t)\bigg|_{t=0}=\frac{1}{2n^2}-\frac{1}{n}log(\theta)+\big(log(\theta)\big)^2. \]
\[ \begin{align} Var\bigg[log\big(X_{(n)}\big)\bigg]&=E\big[log(X_{(n)})^2 \big]- \bigg(E\big[log(X_{(n)}) \big]\bigg)^2\\ &=\frac{1}{2n^2}-\frac{1}{n}log(\theta)+\big(log(\theta)\big)^2- \bigg(-\frac{1}{2n}+log(\theta)\bigg)^2\\ &=\frac{1}{4n^2},\quad \forall\;n\in\mathbb{Z_{\geq 1}}. \end{align} \]
O \(k\)-ésimo momento de \(log\big(X_{(n)}\big)\)
Definição Integral (Gama Incompleta Superior) \[ \Gamma\big(s,x)=\int\limits_{x}^{\infty}t^{s-1}e^{-t}dt, \quad \forall x>0 \quad\wedge \quad s>0. \] Definição em séries de potência (Gama Incompleta Superior) \[ \Gamma\big(s,x)=(s-1)!\;e^{-x}\sum\limits_{k=0}^{s-1}\frac{x^k}{k!}, \quad \forall\; s\in \mathbb{N}\quad \wedge \quad x>0. \] Obs.1): Aqui considerei somente valores reais positivos de \(s\) e \(x\), mas as funções gama incompletas se extendem para outros valores.
Obs.2:) Existe uma relação com a função gama incompleta inferior dada por: \[ \Gamma(x)=\gamma(s,x)+\Gamma(s,x), \quad \forall\;x>0. \]
\[ \begin{align} E\bigg[log\big(X_{(n)}\big)^k\bigg]&=\int\limits_{0}^{\theta}\big(log(x)\big)^k\frac{2n}{\theta^{2n}}x^{2n-1}dx\\ &=\frac{2n}{\theta^{2n}}\int\limits_{0}^{\theta}\big(log(x)\big)^kx^{2n-1}dx\\ &=\frac{(-2n)^{-k}}{\theta^{2n}}\int\limits_{-2nlog(\theta)}^{\infty}y^{k}e^{-y}dy \qquad (Pela\;\;subs. \;\;y=-2nlog(x))\\ &=\frac{(-2n)^{-k}}{\theta^{2n}} \Gamma\big(k+1,-2nlog(\theta)\big), \qquad \forall\; \theta\in(0,1] \quad \wedge \quad k\in \mathbb{Z_{\geq 1}}, \end{align} \]
em que \(\Gamma(a,x)\) é a função gama incompleta superior de parâmetro \(a\) aplicada no ponto \(x>0\).
var.sim2<-function(theta,N,n){
est.umvu<-numeric(length(theta))
for (j in 1:length(theta)){
x_max<-numeric(N)
for (i in 1:N){
u<-runif(n,0,1)
x<-theta*sqrt(u) # amostra gerada (observada)
x_max[i]<-max(x)
}
est.umvu[j]<-var( log(x_max) + 1/(2*n) ) # estimador umvu
}
return(est.umvu)
}
var.teorica<-function(n){1/(4*n^2)}
### Momentos em termos da funcao Gama Incompleta superior
## 1) Implementacao do k-esimo momento usando o pacote zipfR
#install.packages("zipfR") # Usando a funcao Igamma do pacote zipfR
library(zipfR)
Mk<-function(n,k,theta){
((-2*n)^(-k)/(theta^(2*n)))*Igamma(k+1,-2*n*log(theta),lower = FALSE)
}
## 2) Implementacao do k-esimo momento sem o uso de pacote
Igamma.s<-function(s,x){
k<-0:(s-1)
res<-factorial(s-1)*exp(-x)*sum( (x^k)/factorial(k) )
return(res)
}
mk<-function(n,k,theta){((-2*n)^(-k)/(theta^(2*n)))*Igamma.s(k+1,-2*n*log(theta))}## Verificacao de valores
setNames(c(Igamma(2,3,lower=FALSE),Igamma.s(2,3)), c("Gama Inc. Sup. c/ zipfR " , " Gama Inc. Sup s/ pacote"))## Gama Inc. Sup. c/ zipfR Gama Inc. Sup s/ pacote
## 0.1991483 0.1991483setNames(c(Mk(2,3,0.6),mk(2,3,0.6)), c("3º Momento c/ zipfR " , " 3º Momento s/ pacote"))## 3º Momento c/ zipfR 3º Momento s/ pacote
## -0.614313 -0.614313## Testando a funcao momento para alguns casos
M2<-Mk(10,2,0.5) # 1o momento
M1<-Mk(10,1,0.5) # 2o momento
variancia<-M2-M1^2 # usando os momentos
setNames(c(var.teorica(10), variancia, var.sim2(0.5,1000,10)),
c("var.teorica (1/4n^2) ", " Por Momentos ", " var. simulada"))## var.teorica (1/4n^2) Por Momentos var. simulada
## 0.002500000 0.002500000 0.002763692M2<-Mk(100,2,0.7) # 1o momento
M1<-Mk(100,1,0.7) # 2o momento
variancia<-M2-M1^2 # usando os momentos
setNames(c(var.teorica(100), variancia, var.sim2(0.7,1000,100)),
c("var.teorica (1/4n^2) ", " Por Momentos ", " var. simulada"))## var.teorica (1/4n^2) Por Momentos var. simulada
## 2.500000e-05 2.500000e-05 2.323696e-05Questão 08
Dica: Defina-se \(e^{Y}:=1+X\) e encontre a densidade de \(Y\).
Questão 09
Solucão item a):
- Na família exponencial, o EMV é solução da equação \[ E\big[T(\mathbf{X})\big]=T(\mathbf{X}), \quad \mbox{em que T(X) é a estatística suficiente.} \]
Como \(T(\mathbf{X})=\sum\limits_{i=1}^{n}X_i\) (Verifique pelo critério da fatoração!) é estatística suficiente para \(\lambda\),
\[ \begin{align} E\big[\sum\limits_{i=1}^{n}X_i\big]&=\sum\limits_{i=1}^{n}X_i\\ \therefore \quad \sum\limits_{i=1}^{n}E\big[X_i\big]&=\sum\limits_{i=1}^{n}X_i\\ \therefore \quad \frac{n}{\lambda}&=\sum\limits_{i=1}^{n}X_i \quad \because \quad \quad X_1,X_2\ldots,X_n \quad \mbox{são id. distribuídas.}\\ \therefore \quad \widehat{\lambda}_{mv}&=\frac{n}{\sum\limits_{i=1}^{n}X_i}, \quad \mbox{i.e,} \quad \widehat{\lambda}_{mv}=\frac{1}{\bar{X}} \end{align} \]
- Pelo Princípio da Invariância para Estimadores de Máximo-Verossimilhança, \[ \widehat{Var[X_1]}_{mv}=\widehat{\bigg(\frac{1}{\lambda}\bigg)}_{mv}=\frac{1}{\hat{\lambda}_{mv}^2}. \]
ou seja, \[ \widehat{Var[X_1]}_{mv}=\frac{1}{(1/\bar{X})^2}=\bar{X}^2. \]
- Podemos encontrar o estimador UMVU a partir do EMV, uma vez que este último é baseado em estatística suficiente.
Para isso, basta calcular a esperança do EMV e corrigir o viés, caso seja tendencioso, conforme segue:
\[ \begin{align} E\big[\bar{X}^2\big]&=Var[\bar{X}]+\big(E[\bar{X}]\big)^2\\ &=\frac{var[X]}{n}+\big(E[X]\big)\\ &=\frac{1}{n\lambda^2}+\frac{1}{\lambda^2}\\ &=\bigg(1+\frac{1}{n}\bigg)\frac{1}{\lambda^2}=\bigg(\frac{n+1}{n}\bigg)\frac{1}{\lambda^2}. \end{align} \]
Portanto, corrigindo o viés, o estimador UMVU de \(Var[X_1]=\frac{1}{\lambda^2}\) será:
\[ \widehat{Var[X_1]}_{umvu}=\bigg(\frac{n}{n+1}\bigg)\bar{X}^2, \]
em que, conforme visto anteriormente, \(\bar{X}^2\) é o EMV.
Solução do item b)
Defina-se \(Z:=\sum\limits_{i=1}^{n}X_i\). Uma vez que \(X_1,X_2\ldots,X_n\) são i.i.d, tem-se que \(Z\sim Gama(n,\lambda)\).
\[ \begin{align} Var\bigg[\bigg(\frac{n}{n+1}\bigg)\bar{X}^2\bigg]&=\bigg(\frac{n}{n+1}\bigg)^2Var\big[\bar{X}^2\big]\\ &=\frac{n^2}{(n+1)^2}\frac{1}{n^4}Var\bigg[\big(\sum\limits_{i=1}^{n}X_i\big)^2\bigg]\\ &=\frac{1}{[n(n+1)]^2}Var\big[Z^2\big]\\ &=\frac{1}{[n(n+1)]^2}\bigg[E[Z^4]-\big(E[Z^2]\big)^2\bigg] \end{align} \]
Pela função geradora de momentos ou função característica da distribuição \(Gama(n,\lambda)\), \[ E[Z^k]=\frac{\Gamma(k+n)}{\lambda^k\;\Gamma(n)}=\frac{n(n+1)\cdots(n+k-1)}{\lambda^{k}}. \quad \forall\;k\geq1 \quad (\mathbf{Prove!}) \] é valida. Logo,
\[ E[Z^4]=\frac{n(n+1)(n+2)(n+3)}{\lambda^4} \quad \mbox{e}\quad E[Z^2]=\frac{n(n+1)}{\lambda^2}. \]
Desta maneira,
\[ \begin{align} Var\bigg[\bigg(\frac{n}{n+1}\bigg)\bar{X}^2\bigg]&=\frac{1}{\big[n(n+1)\big]^2}\bigg[E[Z^4]-\big(E[Z^2]\big)^2\bigg]\\ &=\frac{1}{\big[n(n+1)\big]^2}\bigg[\frac{n(n+1)(n+2)(n+3)}{\lambda^4} - \bigg(\frac{n(n+1)}{\lambda^2}\bigg)^2\bigg]\\ &=\frac{2(2n+3)}{n(n+1)}\frac{1}{\lambda^4}=\frac{4}{n}\bigg(\frac{n+3/2}{n+1}\bigg) \frac{1}{\lambda^4}\geq \frac{4}{n}\frac{1}{\lambda^4}=LICR \quad \because \quad \frac{n+3/2}{n+1}\geq 1, \quad \forall\;n\geq 1 \end{align} \]
library(ggplot2)
g1<-ggplot() + xlim(1,3) +
geom_function(fun = ~ (2*(2*5 + 3)/(5*(5+1)))*(1/(.x)^4),aes(colour = "Var.UMVU") ) +
geom_function(fun = ~ (4/5)*(1/(.x)^4),aes(colour = "LICR") ) +
labs(title="Variância X LICR", x=bquote(lambda~(n==5)), y="Valores") +
theme(plot.title = element_text(h=0.5))
g2<-ggplot() + xlim(1,3) +
geom_function(fun = ~ (2*(2*20 + 3)/(20*(20+1)))*(1/(.x)^4),aes(colour = "Var.UMVU") ) +
geom_function(fun = ~ (4/20)*(1/(.x)^4),aes(colour = "LICR") ) +
labs(title="Variância X LICR", x=bquote(lambda~(n==20)), y="Valores") +
theme(plot.title = element_text(h=0.5))
library(ggpubr)
ggarrange(g1, g2, ncol=2, nrow=1, common.legend = TRUE, legend="bottom")
library(ggplot2)
g1<-ggplot() + xlim(1,10) +
geom_function(fun = ~ (2*(2*(.x) + 3)/((.x)*((.x)+1)))*(1/0.8^4),aes(colour = "Var.UMVU") ) +
geom_function(fun = ~ (4/(.x))*(1/0.8^4),aes(colour = "LICR") ) +
labs(title="Variância X LICR", x=bquote(n~(lambda==0.8)), y="Valores") +
theme(plot.title = element_text(h=0.5))
g2<-ggplot() + xlim(10,100) +
geom_function(fun = ~ (2*(2*(.x)+3)/((.x)*((.x)+1)))*(1/0.8^4),aes(colour="Var.UMVU") ) +
geom_function(fun = ~ (4/(.x))*(1/0.8^4),aes(colour = "LICR") ) +
labs(title="Variância X LICR", x=bquote(n~(lambda==0.8)), y="Valores") +
theme(plot.title = element_text(h=0.5))
library(ggpubr)
ggarrange(g1, g2, ncol=2, nrow=1, common.legend = TRUE, legend="bottom")
Questão 13
## 1o Modo
## Geracao de amostra
theta <- 2
n <- 10
set.seed(23/03/2021)
x <- rbeta(n,theta,1)
## Verossimilhanca
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,xlab=bquote(theta),type="l",col=4)
## priori
lambda <- 2
r.par <- 2
priori <- dgamma(theta.vec,r.par,rate=lambda)
lines(theta.vec,priori,lty=2)
## Posteriori
posteriori <- dgamma(theta.vec,(n+r.par),rate=(lambda-sum(log(x))))
lines(theta.vec,posteriori,lty=3,col=2)
legend("topright", c("priori", "posteriori", "verossimilhança"),
lty = c(2,3, 1), col = c(1, 2, 4), lwd = c(1.5, 1, 1))
## 2o Modo
set.seed(1234) # fixando semente
n<-10 ; theta<-2 # valores para geracao da amostra
x<-rbeta(n,theta,1) # geracao da amostra (observada)
# verossimilhanca
veross<-function(theta){ sapply(theta,function(theta) prod(dbeta(x,theta,1))) }
# priori
priori <- function(theta, r.par, lambda){dgamma(theta,r.par,rate=lambda)}
# posteriori
posteriori<-function(theta,r.par,lambda){
dgamma(theta,(n+r.par),rate=(lambda-sum(log(x))))}
# Graficos veross, priori e posteriori
library(ggplot2)
ggplot() + xlim(0.01,5)+
geom_function( fun = veross, aes(colour= "veross")) +
geom_function( fun = priori, args =list(r.par=2, lambda=2), aes(colour= "priori")) +
geom_function( fun = posteriori, args=list(r.par=2, lambda=2), aes(colour= "posteriori"))+
labs(title="densidades", x=bquote(theta), y="Valores")+
theme(plot.title = element_text(size=12,hjust=0.5))
Questão 16
priori<-function(alpha,beta,theta){(alpha*beta^alpha)/theta^(alpha+1)}
posteriori<-function(alpha,theta,n){
m<-length(theta)
f<-numeric(m)
for (j in 1:m) {
set.seed(12345)
x.n<-max(runif(n,0,theta[j]))
M<-max(x.n,1)
f[j]<-((n+alpha)/(theta[j]^(n+alpha+1)))*M*ifelse(theta[j]>M,1,0)
}
return(f)
}
veross<-function(theta,n){
m<-length(theta) ;
veros<-numeric(m)
for (j in 1:m) {
set.seed(12345)
x<-runif(n,0,theta[j])
x.n<-max(x)
x.1<-min(x)
veros[j]<-(1/(theta[j]^n))*ifelse(x.1>0 & x.1 < x.n,1,0)*ifelse(theta[j]> x.n,1,0)
}
return(veros)
}
library(ggplot2)
ggplot()+xlim(1.1,2)+
geom_function( fun = priori, args= list(alpha=3,beta=1),aes(colour="Priori"))+
geom_function( fun = veross, args= list(n=5), aes(colour="Verossimilhança"))+
geom_function( fun = posteriori, args = list(alpha=2,n=5), aes(colour="Posteriori"))+
labs(title="densidades", x=bquote(theta), y="Valores")+
theme(plot.title = element_text(size=12,hjust=0.5))
Estimador de Bayes
est.bayes<-function(theta,alpha,n){
m<-length(theta)
est<-numeric(m)
if (all(theta>1)) {
for (j in 1:m){
set.seed(31/03/21)
x<-runif(n,0,theta[j])
x.n<-max(x)
M<-max(x.n,1)
est[j]<-((n+alpha)/(n+alpha-1))*M
}
return(est)
} else {
print("Todo valor de theta deve ser maior que 1 ")
}
}## Estimativas para o parametro theta=3.
n<-c(10,seq(100,1000,100))
bayes.est<-sapply(n, est.bayes,theta=3,alpha=2)
(df<-data.frame(n, bayes.est))## n bayes.est
## 1 10 3.091665
## 2 100 3.005181
## 3 200 2.992868
## 4 300 2.987946
## 5 400 2.995684
## 6 500 2.994197
## 7 600 2.993204
## 8 700 2.992495
## 9 800 3.002109
## 10 900 3.001694
## 11 1000 3.002789Sugestão de Exercício:
Para \(\theta>1\), encontre analiticamente e computacionalmente os erros quadráticos médio:
\[EQM\big(\widehat{\Theta}_{bayes}\big) \quad \mbox{e}\quad EQM\big(\widehat{\Theta}_{mv}\big)\]
E, compare-os.
Questão 18
a<-1; b<-2 #-------------------------- parametros da Beta
delta<-3 ; beta<-5 #------------------ parametros da Gama
n<-1e4 #----------------------------- tamanho da amostra
p<-rbeta(1,a,b) #--------------------- gerando um p (priori beta)
lambda<-rgamma(1,delta,beta) #---------gerando um lambda (priori gama)
x<-numeric(n) ; y<-numeric(n)
for (i in 1:n) {
y[i]<-rpois(1,lambda) # geracao da amostra da dist Poisson
x[i]<-rbinom(1,y[i],p) # geracao da amostra da dist condicional
}
est.p<-(a + sum(x)) / (a + b + sum(y)) # estimativa de Bayes de p
est.lambda<-(delta + sum(y)) / (beta + n) # estimativa de Bayes de lambdasetNames(c(p,est.p),c("Parâmetro p "," Est. de Bayes de p"))## Parâmetro p Est. de Bayes de p
## 0.5803571 0.5861607setNames(c(lambda,est.lambda),c("Parâmetro lambda "," Est. de Bayes de lambda"))## Parâmetro lambda Est. de Bayes de lambda
## 0.6637958 0.6716642Questão 19
\[ \begin{align} \pi(\theta|\mathbf{x})&=\frac{f(\mathbf{x}|\theta)\pi(\theta)}{\int_{\Theta}f(\mathbf{x}|\theta)\pi(\theta)d\theta}\\ &=\frac{\big(\theta c^{\theta}\big)^n\big(\prod\limits_{i=1}^{n}x_i\big)^{-(\theta+1)} \frac{\beta^\delta}{\Gamma(\delta)}\theta^{\delta-1}e^{-\beta\theta}}{\int\limits_{0}^{\infty}\big(\theta c^{\theta}\big)^n\big(\prod\limits_{i=1}^{n}x_i\big)^{-(\theta+1)} \frac{\beta^\delta}{\Gamma(\delta)}\theta^{\delta-1}e^{-\beta\theta}d\theta}\\ &=\frac{\theta^{\delta+n-1}e^{-\theta\big(\beta+\sum\limits_{i=1}^{n}log(x_i/c)\big)}}{\int\limits_{0}^{\infty}\theta^{\delta+n-1}e^{-\theta\big(\beta+\sum\limits_{i=1}^{n}log(x_i/c)\big)}d\theta}\\ &=\frac{\big(\beta+\sum\limits_{i=1}^{n}log(x_i/c)\big)^{\delta+n}}{\Gamma(\delta+n)}\;\theta^{\delta+n-1}\;e^{-\theta\big(\beta+\sum\limits_{i=1}^{n}log(x_i/c)\big)}, \quad \theta>0. \end{align} \] Como \[ F_X(x,\theta)=\int\limits_{-\infty}^{x}f(t)dt=\int\limits_{c}^{x}\theta c^{\theta}t^{-(\theta+1)}dt=1-c^{\theta}x^{-\theta}, \quad \theta>0, \] podemos gerar amostras da distribuição de Pareto do enunciado, a partir da inversa desta função de distribuição acumulada:
\[ F_{X}^{-1}(x,\theta)=c(1-x)^{-1/\theta}, \quad \theta>0. \]
## Parametros iniciais e priori
delta<-2
beta<-3
theta<-rgamma(1,delta,beta)
## geracao de amostra da distribuicao Pareto
n<-1e4 ; c<-2
u<-runif(n,0,1)
x<-c*(1-u)^{-1/theta}
## estimativa de Bayes
theta.hat<-(delta+n) / (beta+sum(log(x/c)))
## Resultado
setNames(c(theta,theta.hat), c("Parâmetro "," Estimativa de Bayes"))## Parâmetro Estimativa de Bayes
## 0.4290748 0.4360105Questão 20
## Parametros inciais e priori
n<-1e4
a<-4; b<-2
theta<-rbeta(1,a,b)
## Amostra de Hardy-Weinberg
y<-sample( c("y1","y2","y3"), size=n, replace=TRUE,
prob=c( theta^2, 2*theta*(1-theta), (1-theta)^2 ))
## Estimativa Bayesiana de theta
theta.hat<- ( 2*sum(y=="y1") + sum(y=="y2") + a ) / (2*n + a + b)## Resultado
setNames(c(theta, theta.hat), c("Parâmetro ", " Estimativa de Bayes"))## Parâmetro Estimativa de Bayes
## 0.8931163 0.8940818## Gerando NPAS da posteriori e calculando a media amostral
amostra.post<-rbeta(n, 2*sum(y=="y1") + sum(y=="y2") + a, sum(y=="y2") + 2*sum(y=="y3") + b)
theta.hat_s<-mean(amostra.post)
setNames(c(theta, theta.hat_s), c("Parâmetro ", " Est. Bayes 2") )## Parâmetro Est. Bayes 2
## 0.8931163 0.8940978