Tarea 7.

Distribución Exponencial

La solución general para la distribución exponencial es

\[l(\theta)=\sum_{j=1}^n(-\ln(\theta)-x_j\theta^{-1})=-n\ln(\theta)-n\overline{x}\theta^{-1},\] \[l'(\theta)=-n\theta^{-1}+n\overline{x}\theta^{-2}=0,\] \[n\theta=n\overline{x},\] \[\hat{\theta}=\overline{x}.\]

Tomando por ejemplo la siguiente data

dsB<-c(27, 82, 115, 126, 155, 161, 243, 294, 340, 384,
457, 680, 855, 877, 974, 1193, 1340, 1884, 2558, 15743)

Se observa gráficamente la función \(l(x)\) con parámetros \(\overline{x}=1424.4\) y \(n=20\).

library(mosaic)
(x_bar<-mean(dsB))
## [1] 1424.4
n_len<-length(dsB)
llexp <- makeFun(-n_len*log(x) - n_len*x_bar/x ~ x)
plotFun(llexp,xlim = range(0,4000))

Distribución Binomial

Para la distribución binomial la solución general es

\(l=\sum_{k=0}^mn_k\Big[\ln\binom{m}{k}+k\ln q +(m-k)\ln(1-q)\Big]\)

Cuando \(m\) es conocida y fija se llega a

\[\hat{q}=\frac{1}{m}\frac{\sum_{k=0}^{\infty}kn_k}{\sum_{k=0}^{\infty}n_k}\].

Con la siguiente data sintética con \(p=0.5\) y \(m=10\) (fija) se observa la verosimilitud.

sample_size <- 30
p <- 0.5
m <- 10
dat1 <- rbinom(sample_size, prob = p, size = m)
table(dat1)
## dat1
## 1 2 3 4 5 6 7 8 
## 1 1 2 5 5 9 6 1
mean(dat1)/10
## [1] 0.5266667
log_likelihood <- function(par_p, par_N, data)
{
   LL <- sum(dbinom(x = data, prob = par_p, size = par_N, log = T))
   return(LL)
}
p_vec <- seq(0, 1, by = 0.01)
logLik <- sapply(p_vec, FUN = log_likelihood, par_N = 10, data = dat1)
par(las = 1, cex.lab = 1.2)
plot(p_vec, logLik, type = "l", xlab = "p", ylab = "lLbinom")
imax <- which.max(logLik)
p_MLE <- p_vec[imax]
p_MLE
## [1] 0.53
points(p_MLE, max(logLik), pch = 19, col = "red")
abline(v = p_MLE, col = "red")

Distribución Poisson

La solución general para la distribución poisson es

\[l=-\lambda n + \sum_{k=0}^{\infty}kn_k\ln\lambda-\sum_{k=0}^{\infty}n_k\ln k!\],

donde \(n=\sum_{k=0}^{\infty}n_k\) es el tamaño de muestra.

Diferenciando con respecto de \(\lambda\)

\[\frac{dl}{d \lambda}=-n + \sum_{k=0}^{\infty}kn_k\frac{1}{\lambda}.\] Igualando con \(0\), el estimador por máxima verosimilitud es

\[\hat{\theta}=\frac{\sum_{k=0}^{\infty}kn_k}{\sum_{k=0}^{\infty}n_k}=\overline{x}.\]

Considerando una data sintética con \(\lambda=2.5\), se tiene la siguiente visualización de la verosimilitud.

sample_size <- 30
lambda<-2.5
dat2 <- rpois(sample_size, lambda)
table(dat2)
## dat2
##  0  1  2  3  4  5 
##  4  4 10  7  2  3
mean(dat2)
## [1] 2.266667
log_likelihood <- function(par_lamb, data)
{
   LL <- sum(dpois(x = data, par_lamb, log = T))
   return(LL)
}
lamb_vec <- seq(0, 5, by = 0.1)
logLik_pois <- sapply(lamb_vec, FUN = log_likelihood, data = dat2)
par(las = 1, cex.lab = 1.2)
plot(lamb_vec, logLik_pois, type = "l", xlab = "lambda", ylab = "lLpois")
imax_pois <- which.max(logLik_pois)
lamb_MLE <- lamb_vec[imax_pois]
lamb_MLE
## [1] 2.3
points(lamb_MLE, max(logLik_pois), pch = 19, col = "red")
abline(v = lamb_MLE, col = "red")

Distribución Geométrica

Con la distribución geométrica, la solución general es

\[\hat{\theta}=\frac{1}{1+\overline{x}},\] al maximizar respecto a la derivada de la verosimilitud

\[\sum_{j=1}^n\Big(\frac{1}{\hat{\theta}}-\frac{x_j}{1-\hat{\theta}}\Big)=0.\]

Se toma la siguiente data con \(p=0.5\) y se observa la verosimilitud.

sample_size <- 30
p <- 0.5
dat3 <- rgeom(sample_size, prob = p)
table(dat3)
## dat3
##  0  1  2  3  4  5 
## 12 11  3  2  1  1
1/(1+mean(dat3))
## [1] 0.483871
log_likelihood_geom <- function(par_p, data)
{
   LL <- sum(dgeom(x = data, prob = par_p, log = T))
   return(LL)
}
p_vec <- seq(0, 1, by = 0.01)
logLik_geom <- sapply(p_vec, FUN = log_likelihood_geom, data = dat3)
## Warning in dgeom(x = data, prob = par_p, log = T): NaNs produced
par(las = 1, cex.lab = 1.2)
plot(p_vec, logLik_geom, type = "l", xlab = "p", ylab = "lLgeom")
imax <- which.max(logLik_geom)
p_MLE <- p_vec[imax]
p_MLE
## [1] 0.48
points(p_MLE, max(logLik_geom), pch = 19, col = "red")
abline(v = p_MLE, col = "red")