Implementação da 1a Prova de Inferência (2/2020)

Modelo \(Gama(2,\theta)\) (Erlang de ordem 2)

lbs<-c("ggplot2","ggpubr")

req <- substitute(require(x, character.only = TRUE))

sapply(lbs, function(x) eval(req) || {install.packages(x,dependencies = TRUE); eval(req)})

## ggplot2  ggpubr 
##    TRUE    TRUE

## funcao de densidade de probabilidade   
densidade <- function(n) {
 layers<-list(ggplot(), 
             xlim(0, 7),  
             xlab("x"),
             ylab("fdp"),
             lapply(1:n, function(theta) 
                   geom_function(fun = ~ theta^2*(.x)*exp(-theta*(.x)), aes(colour=theta))))
 Reduce(ggplot2:::`+.gg`, layers)
}


## funcao de distribuicao acumulada
dist.acum <- function(n) {
 layers<-list(ggplot(), 
              xlim(0, 7),  
              xlab("x"),
              ylab("FDA"),
              lapply(1:n, function(theta) 
                    geom_function(fun = ~ 1-(1+theta*(.x))*exp(-theta*(.x)),aes(colour=theta))))
  Reduce(ggplot2:::`+.gg`, layers)
}

ggarrange(densidade(4), dist.acum(4), ncol=2, common.legend = TRUE, legend="bottom")

Estimadores

Est.mm<-function(amostra){2/mean(amostra)} #Item 1)

Est.mv<-Est.mm  #Item 2)

Est.til<-function(amostra){((2*length(amostra)-1)/(2*length(amostra)))*Est.mv(amostra)} #Item 6) 

Est.chapeu<-function(amostra){mean(ifelse(amostra>1,1,0))} #Item 8)

N<-1e6 ; 

theta<-3 # parametro

x<-rgamma(N,2,theta) # amostra observada da dist gama(2,3)

setNames(c(theta, Est.mm(x),Est.mv(x),Est.til(x)),
         c("theta "," estimativa.mm "," estimativa.mv "," estimativa.til"))

##          theta   estimativa.mm   estimativa.mv   estimativa.til 
##        3.000000        3.001509        3.001509        3.001507

theta<-3

p<-(1+theta)*exp(-theta)

setNames(c(p, Est.chapeu(x)), c("p=P(X>1) "," estimativa de p=P(X>1)"))

##               p=P(X>1)   estimativa de p=P(X>1) 
##               0.1991483               0.1991410

Esperança

Esp.mm<-function(amostra,theta){(2*length(amostra)*theta)/(2*length(amostra)-1)}

Esp.mv<-Esp.mm

Esp.til<-function(amostra,theta){((2*length(amostra)-1)/(2*length(amostra)))*Esp.mv(amostra,theta)}

Vício (Viés)

Vies.mm<-function(amostra,theta){Esp.mm(amostra,theta)-theta}

Vies.mv<-Vies.mm

Vies.til<-function(amostra,theta){Esp.til(amostra,theta)-theta}

setNames(c(Vies.mm(x,3), Vies.mv(x,3),Vies.til(x,3) ) ,
         c("Vicio do e.mm "," Vicio do e.mv "," Vicio do e.til "))

##   Vicio do e.mm    Vicio do e.mv   Vicio do e.til  
##     1.500001e-06     1.500001e-06     0.000000e+00

Erro Quadrático Médio

EQM.mm<-function(n,theta){((n+1)*theta^2)/((n-1)*(2*n-1))}

EQM.mv<-EQM.mm

EQM.til<-function(n,theta){(1/(2*(n-1)))*theta^2}

EQM.chapeu<-function(n,theta){(1+theta)*exp(-theta)*(1-(1+theta)*exp(-theta))/n}

EQM Teórico do EMV X EQM Simulado do EMV

N<-100 ; theta<-3

media.inv<-replicate(1e6,1/mean(rgamma(N,2,3)))

eqm.sim<-4*(N/(N-1))*var(media.inv)+(2*mean(media.inv)-theta)^2

setNames(c(EQM.mv(100,3),eqm.sim),c("EQM Teórico (EMV) "," EQM Simulado (EMV)"))

##  EQM Teórico (EMV)   EQM Simulado (EMV) 
##          0.04613979          0.04657447

n<-50; theta<-seq(2,10,length.out=1000)
df1<-data.frame(X=theta,Y=EQM.mv(n,theta),Z=EQM.til(n,theta))
tipo<-c("Est.mv"="blue","Est.til"="red")

g1<-ggplot() +
    geom_line(data = df1, aes(x = X, y = Y, color = "Est.mv")) +
    geom_line(data = df1, aes(x = X, y = Z, color = "Est.til")) +
    labs(x = bquote(theta ~"( n=50 )"),
       y = "EQM",
       title="Erro Quadrático Médio",
       color = bquote(hat(theta))) +
    scale_color_manual(values = tipo)


theta<-3 ; n<-seq_len(20)
df2<-data.frame(X=n,Y=EQM.mv(n,theta),Z=EQM.til(n,theta))

g2<-ggplot() +
    geom_line(data = df2, aes(x = X, y = Y, color = "Est.mv")) +
    geom_line(data = df2, aes(x = X, y = Z, color = "Est.til")) +
    labs(x = bquote(n~(theta==3)),
         y = "EQM",
         title="Erro Quadrático Médio",
         color = bquote(hat(theta))) +
    scale_color_manual(values = tipo)

ggarrange(g1, g2, ncol=2, common.legend = TRUE, legend="bottom")