## Instalando e Carregando os pacotes (bibliotecas)
lbs<-c('dplyr','ggplot2','ggpubr','DT','kableExtra','fitdistrplus')

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

sapply(lbs, function(x) eval(req) || {install.packages(x,dependencies = TRUE); eval(req)})
##        dplyr      ggplot2       ggpubr           DT   kableExtra fitdistrplus 
##         TRUE         TRUE         TRUE         TRUE         TRUE         TRUE

Questão 03

# geracao dos dados
set.seed(54321)
theta<-4 ; n<-1e4
y<- -log(1-rbeta(n,1,theta))  

# estimativa do parametro
fit <- fitdistr(y, "exponential") 
# Teste de Ajuste ( res.: p-value > 0.05 => distribuicao exp. nao rejeitada )
ks.test(y, "pexp", fit$estimate) 
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  y
## D = 0.0052014, p-value = 0.9496
## alternative hypothesis: two-sided
# plote do grafico
plotdist(y, histo = TRUE, demp = TRUE)

hist(y, freq = FALSE, breaks = 100, xlim = c(0, quantile(y, 0.99)))
curve(dexp(x, rate = fit$estimate), from = 0, col = "red", add = TRUE)

Questões 04 e 06

## LICR e variancia do estimador umvu 
LICR<-function(theta,n){theta^2/n}
var.umvu<-function(theta,n){theta^2/(n-2)}

## Erro Quadratico Medio dos estimadores umvu e mv.
eqm.umvu<-var.umvu
eqm.mv<-function(theta,n){(theta^2/(n-2))*(n/(n-1))^2 + (theta/(n-1))^2}

## Graficos

# Variancia do est. UMVU e LICR
g1<-ggplot() + 
  xlim(1,4) +
  geom_function(fun = LICR , args=list(n=10), aes(colour="LICR")) +
  geom_function(fun = var.umvu , args=list(n=10), aes(colour="Var. UMVU")) +
  labs(title="Variância X LICR", x=bquote(theta~(n==10)), y="Variância") +
  theme(plot.title = element_text(h=0.5))

g2<-ggplot() + 
  xlim(3,20) +
  geom_function(fun = LICR , args=list(theta=4), aes(colour="LICR")) +
  geom_function(fun = var.umvu , args=list(theta=4), aes(colour="Var. UMVU")) +
  labs(title="Variância X LICR", x=bquote(n~(theta==4)), y="Variância") +
  theme(plot.title = element_text(h=0.5))

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

# Erro Quadratico Medio dos estimadores umvu e mv
g3<-ggplot() + 
  xlim(1,20) +
  geom_function(fun = eqm.umvu , args=list(n=10), aes(colour="umvu")) +
  geom_function(fun = eqm.mv , args=list(n=10), aes(colour="mv")) +
  labs(title="Erro Quadrático Médio", x=bquote(theta~(n==10)), y="eqm") +
  theme(plot.title = element_text(h=0.5))

g4<-ggplot() + 
  xlim(3,20) +
  geom_function(fun = eqm.umvu , args=list(theta=4), aes(colour="umvu")) +
  geom_function(fun = eqm.mv , args=list(theta=4), aes(colour="mv")) +
  labs(title="Erro Quadrático Médio", x=bquote(n~(theta==4)), y="eqm") +
  theme(plot.title = element_text(h=0.5))

ggarrange(g3, g4, ncol=2, nrow=1, common.legend = TRUE, legend="bottom")

Questão 07

set.seed(1234)        # fixando semente

n<-10 ; theta<-2      # valores para geracao da amostra

x<-rbeta(n,1,theta)   # geracao da amostra (observada)


## distribuicao a priori
priori<-function(theta,r.par,lambda){dgamma(theta,r.par,lambda)}

## verossimilhanca
veross<-function(theta){ sapply(theta,function(theta) prod(dbeta(x,1,theta))) }

## distribuicao a posteriori
posteriori<-function(theta,r.par,lambda){
                       dgamma(theta,(n+r.par),rate=(lambda+sum(-log(1-x))))}

##  Grafico com  veross, priori e posteriori
ggplot() + xlim(0.01,5)+
      geom_function( fun = veross, aes(colour= "veross")) +
      geom_function( fun = priori, args =list(r.par=4, lambda=2), aes(colour= "priori")) +
      geom_function( fun = posteriori, args=list(r.par=4, lambda=2), aes(colour= "posteriori"))+
      labs(title="densidades", x=bquote(theta), y="Valores")+
      theme(plot.title = element_text(size=12,hjust=0.5))

Questões 05, 08 e 09

## Hiperparametros
r.par<-4 ; lambda<-2  

## Geracao de Amostras
n<-1e4
theta<-rgamma(10,r.par,lambda)
m<-length(theta)
Matrix_sample<-matrix(0,nrow=n,ncol=length(theta))

for (j in 1:m){
  set.seed(12345)
  Matrix_sample[,j]<-rbeta(n,1,theta[j])
}

## Estimadores
est.mv<-function(x){n/sum(- (log(1-x)))}

est.umvu<-function(x){(n-1)/sum(- (log(1-x)))}

est.bayes<-function(x){(n+r.par) / (lambda+sum(-log(1-x))) }

bayes<-apply(Matrix_sample,2,est.bayes)
mv<-apply(Matrix_sample,2,est.mv)
umvu<-apply(Matrix_sample,2,est.umvu)

# Tabela com as estimativas
df<-data.frame(theta,bayes,mv,umvu)
colnames(df)<-c("Parâmetro theta","Est. Bayes", "Est. MV", "Est. UMVU")

kable(df) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
Parâmetro theta Est. Bayes Est. MV Est. UMVU
1.810815 1.817107 1.817041 1.816859
2.762603 2.787531 2.787971 2.787692
1.648358 1.653900 1.653785 1.653620
1.767176 1.776307 1.776227 1.776049
2.471866 2.490423 2.490667 2.490418
1.197098 1.192259 1.192066 1.191947
1.877676 1.885737 1.885694 1.885505
1.321102 1.315399 1.315219 1.315088
1.754687 1.760353 1.760268 1.760092
1.714977 1.725569 1.725475 1.725302