Baseado na página 6 da apostila Verossimilhança e Maxima Verossimiulhança do professor João Luís F. Batista: http://cmq.esalq.usp.br/BIE5781/lib/exe/fetch.php?media=leituras:verossim.pdf

Legenda

Legenda

Link Google

Monitores vs Ocupação

Está abordagem observa o resultado da taxa de ocupação em função do número de orientadores de AAE.
Dataser.

x <- c(11.15, 9.50, 6.65, 14.73, 9.39, 6.23, 23.63)
y <- c(31, 33, 29, 29, 15, 13, 14)

Verificamos a correlação entre os dados.

cor(x, y)
## [1] -0.2243521

Verificação de boxplots.

par(mfrow=c(1, 2))  # divide o grafico em duas colunas

boxplot(x, main="Autuações", sub=paste("Outliers em monitores: ", paste(boxplot.stats(x)$out, collapse=" ")))  # box plot monitores

boxplot(y, main="Ocupação", sub=paste("Outliers na ocupação: ", paste(boxplot.stats(y)$out, collapse=" ")))  # box plot ocupação

Verificação d densidade.

library(e1071)

par(mfrow=c(1, 2))  # dividimos  a área grafica em duas colunas

# Plotamos a densidade do numero de monitores

plot(density(x), main="Densidade: Autuações", ylab="Frequencia", sub=paste("Assimetria:", round(e1071::skewness(x), 2)))

polygon(density(x), col="red")

# Plotamos a densidade da taxa de Ocupação

plot(density(y), main="Densidade: Ocupação", ylab="Frequencia", sub=paste("Assimetria:", round(e1071::skewness(y), 2)))

polygon(density(y), col="red")

Plotamos os dados.

plot(y,x,col = "blue",main = "Relação Autuações vs Oucpação na AAE",
     abline(lm(x~y)),cex = 1.3,pch = 16,
     xlab = "Taxa de Ocupação %",ylab = "Autuações x 10^3")

Verificamos as métricas.

relation <- lm(x~y)
print(summary(relation))
## 
## Call:
## lm(formula = x ~ y)
## 
## Residuals:
##      1      2      3      4      5      6      7 
##  0.682 -0.666 -4.120  3.960 -3.494 -6.956 10.595 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  15.1495     7.2893   2.078   0.0923 .
## y            -0.1510     0.2934  -0.515   0.6286  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.425 on 5 degrees of freedom
## Multiple R-squared:  0.05033,    Adjusted R-squared:  -0.1396 
## F-statistic: 0.265 on 1 and 5 DF,  p-value: 0.6286
 layout(matrix(1:4,2,2))
 plot(relation)

Dispersão.

library(tidyverse)
## Warning: package 'tidyverse' was built under R version 3.4.4
## -- Attaching packages -------------------------------------------------------------- tidyverse 1.2.1 --
## v ggplot2 2.2.1     v purrr   0.2.4
## v tibble  1.4.2     v dplyr   0.7.4
## v tidyr   0.8.0     v stringr 1.3.0
## v readr   1.1.1     v forcats 0.3.0
## Warning: package 'tibble' was built under R version 3.4.4
## Warning: package 'tidyr' was built under R version 3.4.4
## Warning: package 'stringr' was built under R version 3.4.4
## Warning: package 'forcats' was built under R version 3.4.4
## -- Conflicts ----------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(broom)
theme_set(theme_classic())
model.diag.metrics <- augment(relation)
head(model.diag.metrics)
ggplot(model.diag.metrics, aes(x, y)) +
  geom_point() +
  stat_smooth(method = lm, se = FALSE) +
  geom_segment(aes(xend = x, yend = .fitted), color = "red", size = 0.3)

Simulação dos dados

Criamos uma amostra de n=10

dados <- c(24, 27, 23, 28, 26, 24, 17, 23, 24, 27)

Verossimilhança das observações

Obtemos o valor da função da verossimilhança de cada x para (mu= 16)

mu = 16
a<- as.data.frame( sapply(dados, dpois, mu), row.names = NULL, optional = FALSE); 
names(a) <- "a"
head(a)

Obtemos o valor da função da verossimilhança de cada x para (mu= 35)

mu = 35
b<- as.data.frame( sapply(dados, dpois, mu), row.names = NULL, optional = FALSE); 
names(b) <- "b"
head(b)

Obtemos o produtório da verossimilhança das obserações

proda <- prod(a$a)
prodb <- prod(b$b)
prodab <- as.data.frame(cbind(proda, prodb))
prodab

Dataset

Criamos um data frame com os resultados obtidos para mu 16 e 35

ab <- as.data.frame(cbind(a$a, b$b))
names(ab) <- c("mu = 16", "mu = 35")
ab
nplantas <- c(24, 27, 23, 28, 26, 24, 17, 23, 24, 27)
graf1 = as.data.frame(cbind(a, b, nplantas))
graf1 = unique(graf1)
graf1 = graf1[order(graf1$nplantas),c(1,2,3)]
graf1
plot(graf1$nplantas, graf1$a, type = "b", frame = FALSE, pch = 19, col = "red", xlab = "Nº Plantas", ylab = "P(x) para mu= 16")

plot(graf1$nplantas, graf1$b, type = "b", frame = FALSE, pch = 19, col = "red", xlab = "Nº Plantas", ylab = "P(x) para mu= 35")

#http://www.leg.ufpr.br/~paulojus/embrapa/Rembrapa/Rembrapase16.html
lik.pois <- function(lambda, dados) { 
       loglik <- function(l, dados) { 
           sum(dpois(dados, lambda = l, log = TRUE)) 
       } 
       sapply(lambda, loglik, dados = dados) 
}
lambda.vals <- seq(16, 28, l = 101) 
loglik <- sapply(lambda.vals, lik.pois, dados = nplantas)
plot(lambda.vals, loglik, ty = "l")

plot(lambda.vals, loglik, type = "l", xlab = expression(lambda), 
       ylab = expression(l(lambda)))

 curve(lik.pois(x, dados = nplantas), 16, 28, xlab = expression(lambda), ylab = expression(l(lambda)))

dev.pois <- function(lambda, dados) { 
       lambda.est <- mean(dados) 
       lik.lambda.est <- lik.pois(lambda.est, dados = dados) 
       lik.lambda <- lik.pois(lambda, dados = dados) 
       return(-2 * (lik.lambda - lik.lambda.est)) 
} 
curve(dev.pois(x, dados = nplantas), 16, 28, xlab = expression(lambda), 
       ylab = expression(D(lambda)))

 curve(dev.pois(x, dados = nplantas), 14, 30, xlab = expression(lambda), 
       ylab = expression(l(lambda)))

Sigma = x e portanto pode ser obtido diretamente:

lambda.est <- mean(nplantas) 

lambda.est
## [1] 24.3
curve(dev.pois(x, dados = nplantas), 14, 30, xlab = expression(lambda), 
     ylab = expression(l(lambda))) 

L.95 <- qchisq(0.95, df = 1) 

abline(h = L.95)

lim.fc <- function(lambda) dev.pois(lambda, dados = nplantas) - L.95 

ic2.lambda <- c(inf = uniroot(lim.fc, c(0, lambda.est))$root, sup = uniroot(lim.fc, 
     c(lambda.est, max(nplantas)))$root) 
ic2.lambda
##      inf      sup 
## 21.37141 27.48465
arrows(ic2.lambda, L.95, ic2.lambda, 0, len = 0.1) 
 text(ic2.lambda, 0, round(ic2.lambda, dig = 2), pos = 1, cex = 0.8, offset = 0.1)

Acrescentamos o produtório das amostras ao final do dataset

a <- c(a$a, proda)
b <- c(b$b, prodb)
col1 = c("Observação 01", "Observação 02", "Observação 03", "Observação 04", "Observação 05", "Observação 06", "Observação 07", "Observação 08", "Observação 09", "Observação 10", "Produtório")
final <- as.data.frame(cbind(col1, a, b))
names(final) <- c("Itens", "mu = 16", "mu = 35")
final
names(prodab) <- c("mu = 16", "mu = 35")
prodab
ifelse((proda > prodb), "mu 16 é mais verossímil",  "mu 35 é mais verossímil")
## [1] "mu 35 é mais verossímil"

Log de Verossimilhança Negativa

datalog = as.data.frame(cbind(-log(proda), -log(prodb)))
datalog
ifelse((datalog$V1 < datalog$V2), "mu 16 é mais verossímil",  "mu 35 é mais verossímil")
## [1] "mu 35 é mais verossímil"

A Dúvida

O produtório da verossimilhança das observações obtidas foi “1.570865e-20” e “1.938949e-20” para mu 16 e 35 respectivamente; no entando, os resultados da apostila são “2.2574e-22” e “2.2500e-22”.
Qual o erro na tentativa de obtenção do produtório?