Implementação Computacional de Inferência 
Prof. Dr. Antônio Eduardo Gomes1
Aluno Adolfo Manoel Dias da Silva2
26/01/2021
Lista 1
Exercício 01
Item c)
est.p<-function(n,p){
X<-rgeom(n,p) # Geracao da amostra
est1<-1/(1+mean(X))
est2<-2*mean(X) / ( ((n-1)/n)*var(X) + mean(X)*(1+mean(X)) )
est3<-(-1+sqrt( 1+4*((n-1)/n)*var(X) ) ) / ( (2*(n-1)/n)*var(X) )
return(t(list(Estimador_1=est1,Estimador_2=est2, Estimador_3=est3)))
}
est.p(1e6,0.3)## Estimador_1 Estimador_2 Estimador_3
## [1,] 0.2996118 0.3001845 0.3000839Item j)
lambda <- 5
n <- 10 # tamanho da amostra
na <- 100 # número de amostras
mat.dados <- matrix(rexp(n*na,lambda),nrow=na)
no.est <- 11
mat.est <- matrix(0,na,no.est)
mat.est[,no.est] <- 1/apply(mat.dados,1,sd)
for (k in 1:(no.est-1)) {
mat.est[,k] <- (n*factorial(k)/apply((mat.dados^k),1,sum))^(1/k) }
vicios <- apply(mat.est,2,mean)-lambda
variancias <- apply(mat.est,2,var)
par(mfrow=c(1,2))
plot(vicios)
plot(variancias)
Exercício 02
Item h) 1o Modo
Est.Gama<-function(na,alpha,beta){
x<-rgamma(na,alpha,beta) ## Gera a amostra de tamanho na
#### Estimativa de MV para alpha
f<-function(par){-digamma(par)+log(par)+mean(log(x))-log(mean(x))}
alpha.est<-uniroot(f, lower = 0.001, upper = 1e3)$root
#### Estimativa de MV para beta
beta.est<-alpha.est/mean(x)
return(t(list(Estimativa_alpha=alpha.est,Estimativa_beta=beta.est)))
}
Est.Gama(1e4,2,3)## Estimativa_alpha Estimativa_beta
## [1,] 2.029661 3.079714Item h) 2o Modo
logveross <- function(par){
-(par[2]^(n*par[1]))/((gamma(par[1]))^n)*
((prod(x))^(par[1]-1))*exp(-par[2]*sum(x))
}
n <- 30
alfa <- 2
beta <- 3
x <- rgamma(n,alfa,beta)
inicio <- c(1.8,2.8)
emv <- optim(inicio,logveross)
emv$par## [1] 1.905095 2.317795Exercício 03
########## Estudo Computacional Via Simulacao de Monte Carlo ##########
Estimador<-function(m,N,p,I){
#( m: tamanho da amostra ; N:qtd de repeticoes )
#( p: probabilidade ; I: indicador do estimador: 1 ou 2)
R.est<-NULL ; set.seed(12345)
if (I==1) {
for(i in 1:N){
R.est<-c(R.est,sum(rbinom(m,1,p))/(m-sum(rbinom(m,1,p))))
}
} else if (I==2){
for(i in 1:N){
R.est<-c(R.est,1/mean(rgeom(m,p)))
}
} else {
stop("I deve ser 1 ou 2")
}
#### Informacoes do Estimador ######
R2.est<-(R.est)^2
var<-mean(R2.est)-(mean(R.est))^2 # Variancia
vies<-mean(R.est)-p # Vies
eqm<-var + vies^2 # Erro Quad. Medio
return(data.frame(var,vies,eqm))
}
Exercício 04
Sigma.Est<-function(n1,n2,mu1,mu2,sigma){
X<-rnorm(n1,mean=mu1,sd=sigma)
Y<-rnorm(n2,mean=mu2,sd=sigma)
est<- sqrt((sum( (X-mean(X))^2 ) + sum( (Y-mean(Y))^2 ))/(n1+n2) )
return(est)
}
# Exemplo: Considere duas amostras de tamanhos 1000 e 800
# de populacoes normais independentes N_1(2,6) e N_2(3,6).
Sigma.Est(1000,800,2,3,6)## [1] 6.025224Exercício 05
Estimadores<-function(n,r,sigma){
# n: tamanho da amostra
# r: media da normal
# sigma: desvio padrao
A<-pi*r^2 # funcao parametrica
X<-rnorm(n,mean=r,sd=sigma)
X.bar<-mean(X)
S2<-sum((X-X.bar)^2)/(n-1)
A.til<-pi*(X.bar^2-S2/n) # Estimador Não Viesado de A
A.chapeu<-pi*X.bar^2 # Estimador Maxim-Veros de A
return(setNames(c(A,A.til,A.chapeu), c('A','A.til','A.chapeu')))
}
Estimadores(10^4,2,3)## A A.til A.chapeu
## 12.56637 12.18636 12.18915Exercício 07
Ex07<-function(n,theta){
# Geracao de amostra atraves da FDA
u<-runif(n,min=0,max=1)
x<-theta^2/(theta-u)
# Estimadores
Est1<-min(x)
Est2<-((n-1)/n)*Est1
# Erro Quadratico Medio
EQM1<-2*theta^2/( (n-2)*(n-1) )
EQM2<-theta^2/( n*(n-2) )
return(list("Est. de Máximo Veross."=c(Est=Est1,EQM=EQM1),
"Est. Não Viesado"=c(Est=Est2,EQM=EQM2)))
}
Ex07(100,2)## $`Est. de Máximo Veross.`
## Est EQM
## 2.0184475214 0.0008245723
##
## $`Est. Não Viesado`
## Est EQM
## 1.9982630462 0.0004081633
Exercício 08
## Geracao de amostra
gen_amostra<-function(n,theta){runif(n,min=0,max=theta)}
## Estimadores
Estimadores<-function(amostra){
n<-length(amostra)
Est.mm1<-2*mean(amostra)
Est.mm2<-sqrt(12*((n-1)/n)*var(amostra))
Est.mv<-max(amostra)
Est.a_mv<-((n+2)/(n+1))*Est.mv
Est.T<-min(amostra)+max(amostra)
return(setNames(c(Est.mm1,Est.mm2,Est.mv, Est.a_mv,Est.T),
c("Est. MM 1 "," Est. MM 2 "," Est. MV ", " Est. a_MV "," Est. T")))
}
## Erro Quadratico Medio
EQMs<-function(n,theta){
eqm.mm<-theta^2/(3*n) # Atenção: Para calcular EQM do Est.mm2,
# usar resultado do exercicio 26
eqm.mv<-(2*theta^2)/( (n+1)*(n+2) )
eqm.a_mv<-(theta/(n+1))^2
eqm.T<-eqm.mv
return(setNames(c(eqm.mm,eqm.mv,eqm.a_mv,eqm.T),
c("EQM MM 1 "," EQM MV "," EQM a_MV ", " EQM T")))
}
## Exemplo
setNames(theta<-6,"Verdadeiro valor do parâmetro")## Verdadeiro valor do parâmetro
## 6amostra<-gen_amostra(1e8,6)
Estimadores(amostra)## Est. MM 1 Est. MM 2 Est. MV Est. a_MV Est. T
## 5.999745 5.999820 6.000000 6.000000 6.000000EQMs(1e8,theta)## EQM MM 1 EQM MV EQM a_MV EQM T
## 1.2e-07 7.2e-15 3.6e-15 7.2e-15Exercício 10
Item e)
theta <- 3 # verdadeiro valor do parametro na populacao
n <- 10 # tamanho de cada amostra
na <- 1000 # numero de amostras geradas
dados <- matrix(rbeta(n*na,theta,1),nrow=na)
# calculo do EMV para cada amostra gerada
est.mv <- -n/apply(log(dados),1,sum)
# calculo do estimador obtido pelo metodo dos momentos para cada amostra gerada
est.mm <- apply(dados,1,mean)/(1-apply(dados,1,mean))
cbind(mean(est.mv),mean(est.mm))## [,1] [,2]
## [1,] 3.329828 3.273824cbind(sd(est.mv),sd(est.mm))## [,1] [,2]
## [1,] 1.213784 1.216836Exercício 11
## Geracao da amostra
gen_amostra<-function(n,m,p){ rbinom(n,m,p) }
## Estimadores: MV e MM
Estimador<-function(amostra,m){setNames( mean(amostra)/m, "EMV e EMM")}
## Erro Quadratico Medio
EQM<-function(p,n,m){setNames(p*(1-p)/(m*n), "EQM")}
## Exemplo
setNames(p<-0.8,"Verdadeiro Valor do Parametro")## Verdadeiro Valor do Parametro
## 0.8m<-20
amostra<-gen_amostra(1e3,m,p)
Estimador(amostra,m)## EMV e EMM
## 0.80195Exercício 13
## Geracao da amostra de tempos de vida dos componentes
gen_amostra<-function(n,lambda){rexp(n,rate=lambda)}
## Estimador do Tempo de Vida Mediano
Estimador<-function(amostra){setNames(mean(amostra)*log(2),"Est. Tempo Mediano")}
## Erro Quadratico Medio do Estimador Tempo de Vida Mediano
eqm<-function(n,lambda){ e<-(1/n)*(log(2)/lambda)^2 ; e}
## Exemplo
amostra<-gen_amostra(100,6)
Estimador(amostra)## Est. Tempo Mediano
## 0.1204163par(mfrow=c(2,2)) ; n<-1:100
for (i in c(0.02,0.04,0.06,0.08)) {
plot(n,sapply(n,eqm,lambda=i),type="l",xlab="n",ylab="EQM",
main=bquote(paste("Erro Quadrático Médio (",~lambda==.(i),")")))
}
Exercício 14
## Geracao de Amostra
gen_amostra<-function(n,theta) rnorm(n,mean=theta,sd=1)
## Estimadores
Estimadores<-function(amostra){
Env<-mean(amostra>0) # estimador nao viesado
Emv<-pnorm(mean(amostra)) # estimador de max. veross.
return(setNames(c(Env,Emv),c("Est. Não Viesado "," Est. Max. Veross.")))
}
## Erro Quadratico Medio
EQM<-function(amostra){setNames(mean(amostra>0)*(1-mean(amostra>0))/length(n),"EQM de Ybar") }
## Exemplo
theta<-2
setNames(pnorm(0,mean=theta,sd=1,lower.tail=FALSE),"Prob (X > 0)")## Prob (X > 0)
## 0.9772499 amostra<-gen_amostra(1e6,theta)
Estimadores(amostra)## Est. Não Viesado Est. Max. Veross.
## 0.9770490 0.9771209 EQM(amostra)## EQM de Ybar
## 0.0002242425Exercício 15
Item e)
lambda <- 3 # verdadeiro valor do parametro na populacao
n <- 10 # tamanho de cada amostra
na <- 1000 # numero de amostras geradas
dados <- matrix(rpois(n*na,lambda),nrow=na)
# calculo do estimador obtido via E(X) para cada amostra gerada
est1 <- apply(dados,1,mean)
# calculo do estimador obtido via Var(X) para cada amostra gerada
est2 <- apply(dados,1,var)
cbind(mean(est1),mean(est2))## [,1] [,2]
## [1,] 3.0201 3.109656cbind(sd(est1),sd(est2))## [,1] [,2]
## [1,] 0.5518247 1.601227Exercício 16
### Geracao da Amostra a partir da fdp dada
Gen_Amostra<-function(n,theta){
u<-runif(n)
x<-theta*sqrt(u)
return(x)
}
### Estimativas de Maximo Verossimilhanca e por Momentos
Estimador<-function(amostra) {
Emv<-max(amostra)
Env<-(3/2)*mean(amostra)
return(setNames(c(Emv,Env),c("Estimador de MV"," Estimador MM")))
}
### Item d): Verificar qual é o menor EQM
Razao<-function(n) (18*(1+2/n)/((n+1)*(2+1/n)^2))
n<-seq_len(50)
plot(n,sapply(n,Razao),ylim=c(0,2),type="l",ylab="R")
title(main=bquote(R(n)==frac(EQM(hat(theta)[MV]),EQM(hat(theta)[MM]))))
abline(h=1,col="red") ; abline(v=c(4,5),col=c("green","blue"))
axis(1,c(seq_len(10),seq(10,50,5)))
# Exemplo
Amostra<-Gen_Amostra(1000,2)
setNames(min(which(Razao(1:50)<=1)),"tamanho mínimo")## tamanho mínimo
## 5Estimador(Amostra)## Estimador de MV Estimador MM
## 1.998861 2.011830Exercício 17
## Geracao da amostra da fdp dada
gen_amostra<-function(n,theta){
if (theta<=1) {
warning("Theta deve ser maior do que 1 para Estimação!")
} else {
u<-runif(n,0,1)
x<-(1/(1-u))^(1/theta)-1
}
return(x)
}
## Estimadores
Estimadores<-function(amostra){
est.mm<-(mean(amostra)+1)/mean(amostra)
est.mv<- (mean(log(1+amostra)))^{-1}
inv_est.mv<-1/est.mv
return(setNames(c(est.mm,est.mv,inv_est.mv),c("Est. MM"," Est. MV"," Est. MV Inv.")))
}
## Exemplo
amostra<-gen_amostra(1e6,2)
Estimadores(amostra)## Est. MM Est. MV Est. MV Inv.
## 1.9953479 1.9962489 0.5009395Exercício 19
## Gerando amostra da distribuicao com fdp dada.
Gen_Amostra<-function(n,theta,lambda){
u<-runif(n)
x<-(log(1-u)/lambda) + theta
return(x)
}
## Estimador Maximo Verossimilhança
# Em construcaoExercício 21
## Geracao de Amostra de Hardy-Weinberg
gen_amostra<-function(n,theta){
if (theta>=1 || theta <= 0) {
stop("Theta deve ser maior que 0 e menor que 1")
}else{
sample(c('T1',"T2","T3"),size=n,
prob=c(theta^2,2*theta*(1-theta),(1-theta)^2),replace = TRUE)
}
}
## Estimadores de Maximo Verossimilhança
Estimador<-function(amostra){
est.theta<-(2*sum(amostra=="T1")+sum(amostra=="T2"))/(2*length(amostra))
est.func_theta<-est.theta/(1-est.theta)
return(setNames(c(est.theta,est.func_theta),c("Est. theta "," Est. Funcao de theta")))
}
## Exemplo
theta<-0.7
setNames(c(theta,theta/(1-theta)),c("Verdadeiro theta "," Verdadeiro theta/(1-theta)"))## Verdadeiro theta Verdadeiro theta/(1-theta)
## 0.700000 2.333333amostra<-gen_amostra(1e6,theta)
Estimador(amostra)## Est. theta Est. Funcao de theta
## 0.6998255 2.3313956Exercício 22
n<-1000
EQM<-function(a){2*(1-3*a+3*a^2)/((n+2*(n+1)^2))}
a<-seq(0,1,length.out=100)
plot(a,EQM(a),xlab= "a (n=1000)", ylab="EQM",main="Erro Quadrático Médio",type="l")
abline(v=0.5,col="red") ; abline(h=min(EQM(a)),col="blue")
#install.packages("ggplot2")
require(ggplot2)## Loading required package: ggplot2
ggplot() +
xlim(0,1) +
geom_function(fun = EQM , color = "blue")
Exercício 23
Item b)
lambda <- 3
cc <- 0.3
n <- 100 # tamanho de cada amostra
na <- 1000 # numero de amostras geradas
dados.x <- matrix(rexp(n*na,lambda),ncol=n)
dados.y <- sign(sign(cc-dados.x)+1)
est.x <- 1/apply(dados.x,1,mean)
est.y <- (1/cc)*log(n/apply(dados.y,1,sum))
cbind(mean(est.x),mean(est.y))## [,1] [,2]
## [1,] 3.026216 1.745178cbind(sd(est.x),sd(est.y))## [,1] [,2]
## [1,] 0.2944726 0.2801646Item c)
lambda <- 3
# vetor com os 99 percentis da distribuicao Exponencial(lambda)
cc <- qexp(seq(0.01,0.99,0.01),lambda)
n <- 100 # tamanho de cada amostra
na <- 1000 # numero de amostras geradas
# vetor contendo as medias de est.y
# para cada um dos valores de C considerados
media <- rep(0,99)
# vetor contendo os desvios-padrao de est.y
# para cada um dos valores de C considerados
dp <- rep(0,99)
for (i in 1:99) {
dados.x <- matrix(rexp(n*na,lambda),ncol=n)
dados.y <- sign(sign(cc[i]-dados.x)+1)
est.y <- (1/cc[i])*log(n/apply(dados.y,1,sum))
media[i] <- mean(est.y)
dp[i] <- sd(est.y)
}
par(mfrow=c(1,2))
plot(cc,media)
abline(h=lambda)
plot(cc,dp)
Exercício 24
Item b)
lambda <- 1
n <- 100 # tamanho de cada amostra
na <- 1000 # numero de amostras geradas
dados.x <- matrix(rpois(n*na,lambda),ncol=n)
dados.y <- 1-sign(sign(dados.x-1)+1)
est.x <- apply(dados.x,1,mean)
est.y <- log(n/apply(dados.y,1,sum))
cbind(mean(est.x),mean(est.y))## [,1] [,2]
## [1,] 1.00202 1.009872cbind(sd(est.x),sd(est.y))## [,1] [,2]
## [1,] 0.1020124 0.1355631Exercício 26
## Resultado do exercicio 26
Ex26<-function(n,sigma){
esp<-(sqrt(2)*gamma(n/2)/(sqrt(n)*gamma((n-1)/2)))*sigma
varian<-(2*sigma^2/(n*gamma( (n-1)/2) )) * (gamma((n+1)/2) - gamma(n/2)^2/gamma((n-1)/2))
eqm<-varian+(esp-sigma)^2
setNames(c(esp,varian,eqm),c("Esperança "," Variância "," EQM"))
}
## Gerando amostra
gen_amostra<-function(n,mu,sigma){ rnorm(n,mean=mu,sd=sigma) }
## Estimador
est.sigma<-function(amostra){sqrt(((n-1)/n)*var(amostra))}
## Exemplo
amostra<-gen_amostra(190,2,5)
est<-est.sigma(amostra)
setNames(est,"Estimativa de sigma")## Estimativa de sigma
## 4.356358Ex26(190,5)## Esperança Variância EQM
## 4.98023281 0.06570222 0.06609296est.variancia<-function(N,n, mu,sigma){
M<-replicate(N,gen_amostra(n,mu,sigma))
Res<-sqrt(((n-1)/n)*apply(M,2,var))
return( ((N-1)/N) *var(Res) )
}
setNames(est.variancia(1e5,190,2,5),"Estimativa via simulação para variancia teorica")## Estimativa via simulação para variancia teorica
## 0.0656576Exemplo 3.1.6
#devtools::install_github("haozhu233/kableExtra")
library(kableExtra)
## funcao para gerar amostra da fdp ou fda.
amostra<-function(theta,n){
u<-runif(n)
x<-(-1+2*sqrt( 1/4-theta*(1/2-theta/4-u) ) ) / theta
return(x)
}
## Amostra gerada no exemplo do livro
amostra2<-c(0.3374,0.9285,0.6802,-0.2139,0.1052,-0.9793,-0.2623,-0.1964,0.5234, -0.0349,-0.6082, 0.7509,0.3424, -0.7010, -0.2605, 0.4077, -0.7435, 0.9862, 0.9704,0.5313)
## valor inicial para thetaj (thetaj0)
xbar<-mean(amostra2)
## Informação de Fisher
I<-function(par) {(1/( 2*par^3)) * (log( (1+par)/(1-par))- 2*par )}
## funcao que verifica a convergencia dada a amostra, N e theta0
simula<-function(N,theta0,x){
Theta<-theta<-theta0
Valores_Theta<-valores_theta<-NULL
Valores_Theta[1]<-valores_theta[1]<-theta0
for (j in 1:N) {
theta<-theta + sum(x/(1+theta*x)) / sum(x*x/(1+theta*x)^2)
valores_theta[j+1]<-theta
Theta<-Theta + sum(x/(1+Theta*x)) / (length(x)*I(Theta))
Valores_Theta[j+1]<-Theta
}
return(data.frame(Interacao=1:(N+1),Convergencia_1=valores_theta,Convergencia_2 =Valores_Theta))
}
dados<-simula(9,xbar,amostra2)
#dados| Interacao | Convergencia_1 | Convergencia_2 |
|---|---|---|
| 1 | 0.1281800 | 0.1281800 |
| 2 | 0.3587367 | 0.3718426 |
| 3 | 0.3511610 | 0.3491560 |
| 4 | 0.3511317 | 0.3513188 |
| 5 | 0.3511317 | 0.3511139 |
| 6 | 0.3511317 | 0.3511333 |
| 7 | 0.3511317 | 0.3511315 |
| 8 | 0.3511317 | 0.3511317 |
| 9 | 0.3511317 | 0.3511317 |
| 10 | 0.3511317 | 0.3511317 |
Universidade de Brasília, xxxxx@yyyy.com↩︎
Universidade de Brasília, adolfomanoel@hotmail.com↩︎