DESIGUALDAD DE CHEBYSHEV

Este resultado es de carácter teórico y proporciona una cota superior para la probabilidad de que una variable aleatoria tome un valor que diste de su media en más de una cierta cantidad \(\epsilon\) arbitraria.

Sea \(X\) una variable aleatoria con media \(\mu\) y varianza finita \(\sigma^2\). Para cualquier número real \(\epsilon>0\) se tiene que,

• \(\textrm{P}(|X-\mu|>\epsilon) \leq \frac{\sigma^2}{\epsilon^2}\)
• \(\textrm{P}(|X-\mu|>\sigma\epsilon) \leq \frac{1}{\epsilon^2}\)
• \(\textrm{P}(|X-\mu|<\sigma\epsilon) \geq 1-\frac{1}{\epsilon^2}\)

Ejemplo). Sea \(\mu=2\), \(\sigma=4\) y \(\epsilon=2\).

\[ \small{ \begin{align*} &\textrm{P}(|X-2|>4\cdot2) \leq \frac{1}{2^2} \\ &\textrm{P}(-8>X-2>8) \leq \frac{1}{4} \\ &\textrm{P}(-6>X>10) \leq \frac{1}{4} \\ &\textrm{P}(X<-6)+ \textrm{P}(X>10) \leq \frac{1}{4} \\ &\textrm{P}(X<-6)+ [1-\textrm{P}(X<10)]\leq \frac{1}{4} \end{align*} } \]

mu      <- 2
sigma   <- 4
epsilon <- 2

1/(epsilon^2)
#[1] 0.25
pnorm(-6, mean = mu, sd = sigma) + (1-pnorm(10,mean = mu, sd = sigma))
#[1] 0.04550026

En general para cualquier \(\mu\), \(\sigma\) y \(\epsilon\), se tiene que,

\[ \small{ \begin{align*} &\textrm{P}(|X-\mu|>\sigma\epsilon) \leq \frac{1}{\epsilon^2} \\ &\textrm{P}(-\sigma\epsilon>X-\mu>\sigma\epsilon) \leq \frac{1}{\epsilon^2} \\ &\textrm{P}(-\sigma\epsilon+\mu>X>\sigma\epsilon+\mu) \leq \frac{1}{\epsilon^2} \\ &\textrm{P}(X<\mu-\sigma\epsilon)+ \textrm{P}(X>\mu+\sigma\epsilon)\leq \frac{1}{\epsilon^2} \end{align*} } \]

El código para comprobar dicho resultado es una generalización del hecho en el ejemplo, y se presenta a continuación.

mu      <- 0
sigma   <- 1
epsilon <- 2
1/(epsilon^2)

A <- mu-(sigma*epsilon)
B <- mu+(sigma*epsilon)
pnorm(A, mean = mu, sd = sigma) + (1-pnorm(B,mean = mu, sd = sigma))

Ahora vamos a repetir dicho proceso para varios valores de épsilon para ver el comportamiento de la probabilidad y su cota superior. Además, realizaremos la diferencia entre la cota y la probabilidad para observar que efectivamente es una cota, y por lo tanto su diferencia es siempre positiva.

mu      <- 2
sigma   <- 4
n       <- 10
epsilon <- seq(from = 1, to = 10, length = n)
Prob    <- NULL

for (i in 1:n) {
  Prob[i] = pnorm(mu-(sigma*epsilon[i]),mu,sigma)+(1-pnorm(mu+(sigma*epsilon[i]),mu,sigma))
}

A <- cbind(Prob, 1/(epsilon^2),(1/epsilon^2)-Prob)
B <- c("Probabilidad","1/(epsilon^2)","Diferencia")
colnames(A) <- B

Observe el comportamiento de la diferencia entre la probabilidad y la cota superior dejando \(\mu\) y \(\sigma\) fijos y \(\epsilon\) cambiando entre \(1\) y \(10\).

par(col.axis="#CD075A", col.lab="#E36515", cex.lab=1.5, fg="#FF1A1A", mar=c(5,5,0.5,0.5), las=1)
plot(epsilon, (1/epsilon^2)-Prob, type = "l", col = "#FFC71A", lwd = 2,
     xlab = bquote(epsilon), ylab = expression((1/epsilon^2)-Prob))

LEY DÉBIL DE LOS GRANDES NÚMEROS

Sea \(\{X_i\}\) una sucesión de variables aleatorias \(\text{iid}\) con \(\mu = E(X_i)\) y \(\sigma^2 = Var(X_i)\) finita (\(<\infty\)), se tiene que, \[\overline{X}_n \overset{p}{\rightarrow} \mu\]

Donde \(\overline{X}_n = \frac{1}{n}\sum_{i=1}^{n}X_i\).

Ejemplo 1).

n <- 6
p <- 1/4
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
Media_Muestral <- NULL
n1 <- seq(from = 1, to = 500, by = 1)
for (i in 1:length(n1)) {
  Muestra <- rbinom(n = n1[i], size = n, prob = p)
  Media_Muestral[i] <- sum(Muestra)/n1[i]
}
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
par(col.axis="#BD1670", col.lab="#1BA1E2", col.main="#FF4D00", fg="#00D146")
plot(n1, Media_Muestral, type = "l", col = "red", las = 1, ylim = c(1,2),
     main = bquote("X ~ Bin("~n==.(n)~","~p==.(p)~")"),
     xlab = "Tamaño de la Muestra (# Valores simulados)",
     ylab = bquote("Media Muestral ("~bar(x)~")"))
abline(h = n*p, col = "#FEBB0E", lwd = 2)

En el ejemplo anterior podemos observar como a medida que \(n\) aumenta, es decir, el tamaño de la muestra crece, la media muestral tiende a converger (en probabilidad) a la línea amarilla, que en este caso señala el valor esperado de \(X\), es decir, \((n\cdot p)\). Es decir, \(6\cdot \frac{1}{4} = 1.5\) para el ejemplo en cuestión.

Ejemplo 2). A continuación, realizaremos lo mismo del ejemplo anterior, solo que ahora con una función simétrica, que en este caso tomaremos una Normal, y una asimétrica, que en este caso será una Poisson, para observar cuál de las dos converge más rápidamente, o si por lo contrario no hay diferencia.

# X ~ Normal
par(mfrow=c(1,2))
set.seed(8080)
Mu <- 0; Varianza <- 1
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
n <- c(seq(from = 10, to = 500, by = 5), 550, 800, 1000)
Media_Muestral <- NULL

for (i in 1:length(n)) {
  Muestra <- rnorm(n = n[i], mean = Mu, sd = sqrt(Varianza))
  Media_Muestral[i] <- sum(Muestra)/n[i]
}
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
par(col.axis="#00D27F", col.lab="#8A0B80", col.main="#74D600", fg="#0392CF")
plot(n, Media_Muestral, type = "l", col = "orange", las = 1,
     main = bquote("X ~ N("~mu==.(Mu)~","~sigma^2==.(Varianza)~")"),
     xlab = "Tamaño de la Muestra (# Valores simulados)",
     ylab = "Media Muestral", ylim = c(-0.4,0.4))
abline(h = Mu, col = "#7044A6", lwd = 2)
# __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ #
# X ~ Poisson
set.seed(8080)
Lambda <- 2
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
n <- c(seq(from = 10, to = 500, by = 5), 550, 800, 1000)
Media_Muestral <- NULL

for (i in 1:length(n)) {
  Muestra <- rpois(n = n[i], lambda = Lambda)
  Media_Muestral[i]  <- sum(Muestra)/n[i]
}
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
plot(n, Media_Muestral, type = "l", col = "orangered", las = 1,
     main = bquote("X ~ Poisson("~lambda==.(Lambda)~")"),
     xlab = "Tamaño de la Muestra (# Valores simulados)",
     ylab = "Media Muestral", ylim = c(1.6,2.4))
abline(h = Lambda, col = "#7044A6", lwd = 2)

DISTRIBUCIÓN LÍMITE DE \(\overline{x}_n\)

\[ \lim_{n\rightarrow \infty}F_n(\overline{x}_n) = \begin{cases} 0 & \textrm{si } \overline{x}_n<0 \\ 1/2 & \textrm{si } \overline{x}_n=0 \\ 1 & \textrm{si } \overline{x}_n>0 \end{cases} \]

FxMedia <- function(n){
  par(mfrow = c(1,2))
  par(col.axis="#BD1670", col.lab="#E32E01", col.main="#FBBE00", fg="#00D146")
  x <- seq(from = -1, to = 1, by = 0.001)
  if (n<1000) {
    D <- dnorm(x, 0, 1/sqrt(n))
    plot(x, D, xlim = c(-0.5,0.5), col = "deepskyblue2", las = 1, lwd = 2,
         type = "l", main = "FUNCIÓN DE DENSIDAD", xlab = "x", ylab = "f(x)")
    curve(pnorm(x, 0, 1/sqrt(n)), xlim = c(-0.5,0.5), col = "darkorange1", las = 1, lwd = 2,
          main = "FUNCIÓN DE DISTRIBUCIÓN", xlab = "x", ylab = "F(x)")
    abline(h = 0.5, col = "firebrick1")
  }
  else{
    D <- dnorm(x, 0, 1/sqrt(n))
    plot(x, D, xlim=c(-1,1), col = "deepskyblue2", las = 1, lwd = 2,
         type = "l", main = "FUNCIÓN DE DENSIDAD", xlab = "x", ylab = "f(x)")
    curve(pnorm(x, 0, 1/sqrt(n)), xlim = c(-1,1), col = "darkorange1", las = 1, lwd = 2,
          main = "FUNCIÓN DE DISTRIBUCIÓN", xlab = "x", ylab = "F(x)")
    abline(h = 0.5, col = "firebrick1")
  }
}
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
FxMedia(100)

FxMedia(5000)

ESTADÍSTICAS DE ORDEN

Vamos a comprobar que, en general,

\[ X_{i} {\color{DarkOrange} {\Rightarrow}} \begin{cases} {\rm I\!E}(X_i) &= \mu \\ \textrm{Var}(X_i) &= \sigma^2 \end{cases} \\ X_{(i)} {\color{red} {\nRightarrow}} \begin{cases} {\rm I\!E}(X_{(i)}) &= \mu \\ \textrm{Var}(X_{(i)}) &= \sigma^2 \end{cases} \]

Para ello vamos a crear la siguiente matriz, en donde cada columna será una variable aleatoria ordenada de menor a mayor (en este caso solamente lo haremos con \(10\)), y cada registro será una muestra (en este caso tomaremos \(1000\)) de cualquier distribución.

En este caso suponga que \(X \overset{d}{=} \mathcal{N}(\mu, \sigma)\), el siguiente código contempla la matriz antes mencionada.

set.seed(54321)
mu    <- 10
sigma <- 1
Tam_Muestra <- 10
Tam_Simu    <- 1000
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# Estadísticas de Orden
EO <- matrix(0, Tam_Simu, Tam_Muestra)
for (i in 1:Tam_Simu) {
  Muestra = rnorm(Tam_Muestra, mu, sigma)
  EO[i,]  = sort(Muestra)
}

Así hemos obtenido una muestra aleatoria de las \(10\) (Tam_Muestra) variables aleatorias ordenadas (estadísticas de orden).

head(EO)

dim(EO)
# [1] 1000   10
EO[1,]
# [1]  8.307758  8.349400  8.904471  9.071956  9.215966  9.591933  9.821099 10.179977 11.395352 12.516046
summary(EO[,1])
#  Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
# 6.137   8.048   8.471   8.430   8.831  10.270

A continuación, graficaremos las densidades de las estadísticas de orden \(X_{(1)}\) y \(X_{(10)}\) (mínimo y el máximo).

par(mfrow=c(1,2))
par(col.axis="#7633F1", col.lab="#830C4C", col.main="#EB015A")
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
hist(EO[, 1], col = ColorF, density = 50, prob = T, las = 1,
     xlab = bquote(X[(1)]), ylab = "Densidad", main = "DENSIDAD DEL MÍNIMO")
lines(density(EO[,1]), col = "#0F81D8", lwd = 2)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
hist(EO[, 10], col = ColorG, density = 50, prob = T, las = 1,
     xlab = bquote(X[(10)]), ylab = "Densidad", main = "DENSIDAD DEL MÁXIMO")
lines(density(EO[,10]), col = "#2BB800", lwd = 2)

Ejemplo 1). Sea \(X_1, X_2, X_3 \textrm{ y } X_4\) una \(\textit{m.a.}\) de \(X \overset{d}{=} \textrm{Beta}(\alpha = 2, \beta = 2)\). Encontrar la probabilidad de que el máximo \(X_{(4)}\) de la muestra sea menor a \(0.9\).

Recuerde que,

• \(F_{X_{(r)}}(x) = \sum_{i=r}^{n}\binom{n}{i}[F_{X}(x)]^i [1-F_{X}(x)]^{n-i}\)
• \(F_{X_{(1)}}(x) = 1-[1-F_{X}(x)]^n\)
• \(F_{X_{(n)}}(x) = [F_{X}(x)]^n\)

\[\begin{align*} P(X_{(4)}<0.9) &= [F_X(x)]^n \\ &= \left( \int_{0}^{0.9}6x(1-x)dx \right)^4 \\ &= (0,972)^4 \\ &= 0,892616 \end{align*}\]

(pbeta(0.9,2,2))^4

## [1] 0.8926168

Ahora vamos a comprobar el resultado anterior mediante simulación.

set.seed(2000)
Tam_Simu    <- 1000
Tam_Muestra <- 4
X4 <- NULL
for ( i in 1:Tam_Simu) {
  muestra = rbeta(Tam_Muestra,2,2)
  X4[i] = max(muestra)
}

par(col.axis="#FF2640", col.lab="#FF0DBF", col.main="#FF3300")
hist(X4, col = ColorB, density = 50, prob = T, las = 1,
     xlab = bquote(X[(4)]), ylab = "Densidad", main = "DENSIDAD DEL MÁXIMO")
lines(density(X4), col = "#1CC513", lwd = 2)

Con ello, procederemos a calcular la probabilidad.

CasoFav <- NULL
for (i in 1:Tam_Simu) {
  CasoFav[i] <- ifelse(test = X4[i]<0.9, yes = 1, no = 0)
}
( Probabilidad <- sum(CasoFav)/Tam_Simu )

## [1] 0.892

Ejemplo 2). Sea \(X_{(n)}\) la n-ésima estadística de orden de una \(\textit{m.a.}\) de \(X \overset{d}{=} \textrm{U}(0, \theta)\). Encontrar la distribución límite de \(X_{(n)}\).

\[ F_X(x)=\left\{\begin{matrix} 0 & x<0 \\ (\frac{x}{\theta}) &0\leq x<\theta \\ 1 & x\geq \theta \end{matrix}\right. \]

\[\begin{align*} F_{X_{(n)}} &= (F_X(x))^{n} \\ &= \left ( \frac{x}{\theta } \right )^{n} \\ \lim_{n\rightarrow \infty}F_{X_{(n)}} &= \lim_{n\rightarrow \infty}\left ( \frac{x}{\theta } \right )^{n} \\ &= \left\{\begin{matrix} 0 & x\leq \theta \\ 1 & x>\theta \end{matrix}\right. \end{align*}\]

Theta <- 1
n <- c(5, 10, 20, 30, 50, 100, 200)

Colors <- c("#BF1029", "#009ABA", "#4FC009", "#BF1369", "#FF6314", "#EEA804", "#1C8F6A")
par(col.axis="#7C00A3", col.lab="#29ABE1", cex.lab=1.5, col.main="#FF4500", cex.main=1.8,
    fg="#D11141", mar=c(5,5,2,0.5), las=1)
for(i in 1:length(n)){
  curve((x/Theta)^n[i], xlim = c(0,Theta),
        add = i>1, type = "l", col = Colors[i], lwd = 2,
        main = "Función de Distribución", xlab = bquote(X[(n)]),
        ylab = bquote((x/theta)^n))
  legend(0, 0.06*i+0.6, legend = bquote(n==.(n[i])),
         col = Colors[i], lty = 1, box.col = "transparent", bg = "transparent",
         text.col = "black", cex = 1.3)
}

DISTRIBUCIONES MUESTRALES

MÉTODO DE LOS MOMENTOS

Por cómo está definido el estimador de momentos al igual que el estimador de máxima verosimilitud, puede que este no sea único. Por ejemplo, tenemos una muestra aleatoria proveniente de una distribución \(\textrm{Pois}(\lambda)\), se sabe que \({\rm I\!E}(X) = \lambda\), entonces un estimador de momentos de \(\lambda\) es, naturalmente, \(\overline{X}\). Pero también se sabe que \(\textrm{Var}(X) = \lambda\), de esta manera, se tiene otro estimador de \(\lambda\) que es \(S_n^2\).

La pregunta natural ahora es ¿cuál de las dos estimaciones es mejor?, es decir, ¿cuál estimador se acerca más al valor verdadero de \(\lambda\)? No podemos responder a esta pregunta directamente, puesto que no conocemos el valor de \(\lambda\). Más adelante veremos conceptos que nos permiten evaluar la calidad de los estimadores, por ahora un simple ejercicio de simulación nos permitirá elegir uno de forma empírica.

Vamos a simular muestras provenientes de una distribución \(\textrm{Pois}(5)\) y \(\textrm{Pois}(15)\) con un tamaño de muestra \(n=5,\cdots ,300\) y para cada muestra simulada, se calculan las dos estimaciones de momentos \(\overline{x}\) y \(s_n^2\), y se observa cuál es más cercano al valor verdadero del parámetro. La línea amarilla horizontal denota el valor verdadero del parámetro \(\lambda\). El comando en \(\texttt{R}\) y los resultados se presentan a continuación.

set.seed(4321)
Mean_Var.Poisson <- function(Lambda) {
  n <- 5:300
  Est.Mean <- NULL
  Est.Var  <- NULL
  for (i in 1:length(n)) {
    Muestra     <- rpois(n[i], Lambda)
    Est.Mean[i] <- mean(Muestra)
    Est.Var[i]  <- (n[i]-1)/n[i]*var(Muestra)
  }
  
  par(col.lab="#A32062", col.main="#FF4D00", fg="#069900", las=1)
  plot(x = n, y = Est.Var, type = "l", col = "#00ABA9",
       xlab = "n", ylab = "Estimación", ylim = c(Lambda-4,Lambda+4),
       main = bquote("X ~ Poisson("~lambda==.(Lambda)~")"))
  lines(n, Est.Mean, col = "#D41243")
  abline(h = Lambda, col = "#FEBB0E", lwd = 2)
  
  legend("bottomright", c("Media","Varianza"), lty = c(1,1),
         col = c("#D41243","#00ABA9"), box.col='transparent',
         bg = 'transparent', text.col = 1, cex = 0.9)
}

par(mfrow=c(2,1))
Mean_Var.Poisson(Lambda = 5); Mean_Var.Poisson(Lambda = 15)

Se puede ver claramente que la media \(\overline{X}\) comparada con \(S_n^2\) estima mejor el parámetro de la distribución, puesto que las estimaciones de \(\overline{X}\) parecen estar más cercanas del valor de \(\lambda\) en ambas gráficas. Más adelante se discutirán métodos formales acerca de la elección entre estimadores, y se verá que el estimador \(\overline{X}\) tiene mejores propiedades que \(S_n^2\).

MÉTODO DE MÁXIMA VEROSIMILITUD

La función de verosimilitud se define como:

\[\mathcal{L}(X|\theta) = \prod_{i=1}^n f(X_i|\theta) = f(X_1|\theta)f(X_2|\theta)\cdots f(X_n|\theta)\]

y la log-verosimilitud se define como:

\[\log\mathcal{L}(X|\theta)= \log\left(\prod_{i=1}^n f(X_i|\theta)\right) = \sum_{i=1}^n \log(f(X_i|\theta)) = \log(f(X_1|\theta))+\cdots+\log(f(X_n|\theta))\]

Ejemplo 1). Sea \(X \overset{d}{=} \textrm{Exp}(\lambda)\) hallar el estimador de máxima verosimilitud de \(\lambda\).

\[\begin{align*} \mathcal{L}(\lambda) &= \prod_{i=1}^{n}\lambda e^{-\lambda x_i} \\ &= \lambda^n e^{-\lambda \sum_{i=1}^{n} x_i} \\ \ln (\mathcal{L}(\lambda)) &= n\ln(\lambda)-\lambda \sum_{i=1}^{n} x_i \\ \frac{\partial \ln (\mathcal{L}(\lambda))}{\partial \lambda} &= \frac{n}{\lambda} - \sum_{i=1}^{n} x_i = 0\\ \frac{n}{\lambda} &= \sum_{i=1}^{n} x_i \\ \frac{n}{\sum_{i=1}^{n} x_i} &= \lambda \\ \frac{1}{\overline{x}} &= \lambda \\ \frac{\partial^2 \ln (\mathcal{L}(\lambda))}{\partial \lambda^2} ) &= -\frac{n}{\lambda^2} \\ \widehat{\lambda} &= \frac{1}{\overline{x}} \end{align*}\]

Observe que la segunda derivada es estrictamente negativa, luego \(\widehat{\lambda}\) es efectivamente un máximo. Con ello tenemos que \(\widehat{\lambda} = \frac{1}{\overline{x}}\), es decir el inverso de la media empírica, lo que es coherente con el hecho que el parámetro \(\lambda\) es el inverso de la esperanza.

Ejemplo 2). Sea \(X \overset{d}{=} \textrm{Geo}(p)\) \(f(x) = p(1-p)^{x-1}\) hallar el estimador de máxima verosimilitud de \(p\).

\[\begin{align*} \mathcal{L}(p) &= \prod_{i=1}^{n} p(1-p)^{x_i-1} \\ &= p^n (1-p)^{\sum_{i=1}^{n} x_i - n} \\ \ln (\mathcal{L}(p)) &= n\ln(p) + (\sum_{i=1}^{n} x_i - n)\ln(1-p) \\ \frac{\partial \ln (\mathcal{L}(p))}{\partial p} &= \frac{n}{p}-\frac{\sum_{i=1}^{n} x_i - n}{1-p} = 0\\ 0 &= \frac{n-np-p\sum_{i=1}^{n} x_i + np}{p(1-p)} \\ 0 &= n-p\sum_{i=1}^{n} x_i \\ p &= \frac{n}{\sum_{i=1}^{n} x_i} = \frac{1}{\overline{x}} \\ \frac{\partial^2 \ln (\mathcal{L}(p))}{\partial p^2} ) &= -\frac{n}{p^2}-\frac{\sum_{i=1}^{n} x_i - n}{(1-p)^2} \\ \widehat{p} &= \frac{1}{\overline{x}} \end{align*}\]

Observe que la segunda derivada es estrictamente negativa, luego \(\widehat{p}\) es efectivamente un máximo. Con ello tenemos que \(\widehat{p} = \frac{1}{\overline{x}}\), es decir el inverso de la media empírica, lo que es coherente con el hecho que el parámetro \(p\) es el inverso de la esperanza.

Ejemplo 3). Sea \(X \overset{d}{=} \mathcal{N}(\mu,\sigma^2)\) hallar el estimador de máxima verosimilitud de \(\mu\) y \(\sigma\).

\[\begin{align*} \mathcal{L}(\mu,\sigma^2) &= \prod_{i=1}^{n} \frac{1}{\sqrt{2\pi\sigma^2}} e^{\frac{-(x_i-\mu)^2}{2\sigma^2}} \\ &= \left ( \frac{1}{\sqrt{2\pi\sigma^2}} \right )^n \text{exp} \left \{ -\frac{1}{2\sigma^2} \sum_{i=1}^{n}(x_i-\mu)^2 \right \} \\ \ln (\mathcal{L}(\mu,\sigma^2)) &= n\ln \left ( \frac{1}{\sqrt{2\pi\sigma^2}} \right ) - \frac{1}{2\sigma^2} \sum_{i=1}^{n}(x_i-\mu)^2 \\ \frac{\partial \ln (\mathcal{L}(\mu,\sigma^2))}{\partial \mu} &= \frac{1}{\sigma^2} \sum_{i=1}^{n}(x_i-\mu) \\ \frac{\partial \ln (\mathcal{L}(\mu,\sigma^2))}{\partial \sigma^2} &= -\frac{n}{2\sigma^2}+\frac{1}{2(\sigma^2)^2}\sum_{i=1}^{n}(x_i-\mu)^2 \\ 0 &= \frac{1}{\sigma^2} \sum_{i=1}^{n}x_i-\mu \\ 0 &= -\frac{n}{2\sigma^2}+\frac{1}{2(\sigma^2)^2}\sum_{i=1}^{n}(x_i-\mu)^2 \end{align*}\]

Los estimadores de máxima verosimilitud de \(\mu\) y \(\sigma\) son respectivamente la media y la varianza empíricas de la muestra, \(\overline{x}\) y \(s_n^2\), tal como era de esperar.

Ejemplo 4). Sea \(X:=\) \(\#\) de accidentes de tránsito en la calle 26 en una semana.

\[X \overset{d}{=} \textrm{P}(\lambda)\]

Si se observó que \(x = 2\), hallar la estimación por máxima verosimilitud de \(\lambda\).

# En una muestra se observó x=2
x <- 2
Lambda <- seq(from = 0.1, to = 5, by = 0.1)
Proba  <- dpois(x, Lambda)

par(col.lab="#FF4D00", cex.lab=1.5, fg="#FFAB19", mar=c(5,5,0.5,0.5), las=1)
plot(Lambda, Proba, type = "l", col = "#FB2D60", lwd = 3,
     xlab = bquote(lambda), ylab = bquote("Probabilidad de "~X==.(x)))
abline(v = x, col = "#24C8A0", lwd = 2, lty = 2)

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
Names <- c("Lambda", "P(X=2)")
A     <- cbind(Lambda, round(Proba,5))
colnames(A) <- Names

# __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ #
# En otra muestra se observó x=3
y <- 3
Proba2 <- dpois(y, Lambda)
# Producto
ProbConjunta <- Proba*Proba2

par(col.lab="#FF4D00", cex.lab=1.5, fg="#FFAB19", mfrow=c(1,2), mar=c(5,5,0.5,0.5), las=1)
plot(Lambda, ProbConjunta, xlab = bquote(lambda), ylab = "P(X=2, X=3)",
     type = "l", col = "#74D600", lwd = 3)
abline(v = (x+y)/2, col = "#DE0042", lwd = 2, lty = 2)
plot(Lambda, log(ProbConjunta), xlab = bquote(lambda), ylab = "P(X=2, X=3)",
     type = "l", col = "#1EBBD7", lwd = 3)
abline(v = (x+y)/2, col = "#DE0042", lwd = 2, lty = 2)

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
max(ProbConjunta)

## [1] 0.05483355

Names <- c("Lambda", "P(X=2,X=3)")
B     <- cbind(Lambda,round(ProbConjunta,5))
colnames(B) <- Names

\(\texttt{R}\) puede usarse para encontrar estimadores de máxima verosimilitud en muchos entornos diversos. Discutiremos solo lo más básico aquí y dejaremos el resto a textos más sofisticados.

Para problemas de estimación de un parámetro, podemos usar la función optimize() para encontrar el EMV (estimador máximo verosímil). Los argumentos son, la función a maximizar (la función de probabilidad), el rango sobre el cual se realizará la optimización y, opcionalmente, cualquier otro argumento que se pasará a la probabilidad si es necesario.

Realicemos como ejemplo, la función de probabilidad conjunta dada anteriormente.

x <- 2:3
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
L <- function(x, Lambda) { prod(dpois(x, lambda = Lambda)) }
optimize(L, interval = c(0,5), x = x, maximum = TRUE)

## $maximum
## [1] 2.5
## 
## $objective
## [1] 0.05483355

Tenga en cuenta que la probabilidad es solo un producto de las funciones de probabilidad de las Poisson. Primero damos algunos datos de muestra (en los valores de datos vectoriales), luego definimos la función de probabilidad L, y finalmente optimizamos L sobre el rango c(0,5), es decir, \(0\leq \lambda \leq5\).

Tenga en cuenta que la función de optimización por defecto minimiza la función L, por lo que tenemos que establecer maximum = TRUE para obtener un EMV. El valor devuelto de $máximo da un valor aproximado del EMV de \(2.5\) y $objective da L evaluado en el EMV que es aproximadamente \(0.05483355\).

Antes comentamos que generalmente es numéricamente más conveniente maximizar la función log-verosimilitud y podemos hacerlo con la misma facilidad en \(\texttt{R}\).

x <- 2:3
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
MinimoLogL <- function(x, Lambda){ -sum(dpois(x, lambda = Lambda)) }
optimize(MinimoLogL, interval = c(0,5), x = x)

## $minimum
## [1] 2.449502
## 
## $objective
## [1] -0.470496

Tenga en cuenta que no necesitábamos maximum = TRUE porque minimizamos la probabilidad logarítmica negativa. La respuesta para el EMV es esencialmente la misma que antes, pero el valor de $objective fue diferente, por supuesto.

COMPARACIÓN ENTRE EL MÉTODO DE LOS MOMENTOS Y EL DE MÁXIMA VEROSIMILITUD

Sea \(X_1, \cdots,X_n\) una muestra aleatoria proveniente de una distribución uniforme continua sobre \((0,\theta)\). En donde para encontrar el estimador de momentos de \(\theta\), se tiene presente que \(\mu = \theta/2\) de donde \(\theta = 2\mu\), entonces se concluye que \(\widehat{\theta}_{mom} = 2\overline{X}\), el cual es diferente al estimador por máxima verosimilitud dado por \(\widehat{\theta}_{MV} = X_{(n)}\)

set.seed(2019)
Unif_Esti <- function(Theta) {
  n <- seq(from = 1, to = 300, by = 1)
  Est_Momen <- NULL
  Est_MV    <- NULL
  
  for (i in 1:length(n)) {
    Muestra       <- runif(n[i], min = 0, max = Theta)
    Est_Momen[i]  <- 2*mean(Muestra)
    Est_MV[i]     <- max(Muestra)
  }
  
  par(col.lab="#A32062", col.main="#FF4D00", fg="#069900", las=1)
  plot(x = n, y = Est_Momen, type = "l", col = "#00ABA9",
       xlab = "n", ylim = c(Theta-1,Theta+1),
       ylab = bquote("Estimación para "~theta),
       main = bquote("X ~ Uniforme( 0,"~.(Theta)~")"))
  
  lines(n, Est_MV, col = "#D41243")
  abline(h = Theta, col = "#FEBB0E", lwd = 2)
  
  legend("bottomright", c("MOMENTOS","MV"), lty = c(1,1),
         col = c("#00ABA9","#D41243"), box.col = 'transparent',
         bg = 'transparent', text.col = 1, cex = 0.9)
}

par(mfrow=c(2,1))
Unif_Esti(Theta = 3); Unif_Esti(Theta = 5)

Podemos ver con el estimador de máxima verosimilitud siempre se obtuvieron valores más cercanos al parámetro sin importar el tamaño muestral, aunque las estimaciones de máxima verosimilitud parecen estar por debajo del \(\theta\) verdadero, situación que no sucede con las estimaciones de momentos.

PROPIEDADES DE LOS ESTIMADORES PUNTUALES

mu    <- 0
sigma <- 2
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
#       Recorrido de la función para graficarla
Min <- mu-3*sigma; Max <- mu+3*sigma
x   <- seq(Min, Max, length = 10000)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
Medias          <- NULL
varianza_Inses  <- NULL
varianza_Ses    <- NULL

Tam_Muestra <- 1000

for (i in 1:Tam_Muestra) {
  Muestra           = rnorm(Tam_Muestra, mu, sigma)
  Medias[i]         = mean(Muestra)
  varianza_Inses[i] = var(Muestra)
  varianza_Ses[i]   = (Tam_Muestra-1)/Tam_Muestra*var(Muestra)
}
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
mean(Medias)

## [1] 0.001375969

mean(varianza_Inses)

## [1] 4.009624

mean(varianza_Ses)

## [1] 4.005614

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
par(col.lab="#A32062", cex.lab=1.3, col.main="#FF4D00", cex.main=1.6,
    fg="#069900", mar=c(4,4,2,0.5), las=1)
hist(Medias, col = rainbow(8), density = 50, prob = T,
     main = bquote("HISTOGRAMA DEL COMPORTAMIENTO DE "~bar(X)),
     ylab = "Densidad", xlab = bquote(bar(X)))
lines(x, dnorm(x,mu,sigma/sqrt(Tam_Muestra)), type = "l", col = 4, lwd = 3)
abline(v = mu, col = 2, lwd = 2, lty = 2)

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
par(mfrow=c(2,1))
hist(varianza_Inses, col = rainbow(8), density = 50, prob = T,
     main = bquote("HISTOGRAMA DEL COMPORTAMIENTO DE "~S^2),
     xlab = bquote(S^2), ylab = "Densidad")
lines(x, dnorm(x, sigma^2 , sqrt(2*sigma^4/(Tam_Muestra-1))),
      type = "l", col = 4, lwd = 3)
abline(v = sigma^2, col = 2, lwd = 2, lty = 2)

hist(varianza_Ses, col = rainbow(8), density = 50, prob = T,
     main = bquote("HISTOGRAMA DEL COMPORTAMIENTO DE "~S[n]^2),
     xlab = bquote(S[n]^2), ylab = "Densidad")
lines(x, dnorm(x, sigma^2 ,sqrt(2*sigma^4/(Tam_Muestra))),
      type = "l", col = 4, lwd = 3)
abline(v = sigma^2, col = 2, lwd = 2, lty = 2)

set.seed(1234)
mu    <- 0
sigma <- 1

Tam_Muestra   <- 10   # Probar con 100 y 1000
Num_Muestras  <- 1000
Varianzas     <- NULL

for (i in 1:Num_Muestras) {
  Muestra      = rnorm(Tam_Muestra, mu, sigma)
  Varianzas[i] = var(Muestra)
}

A <- ((Tam_Muestra-1)*Varianzas)/(sigma^2)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
#       Recorrido de la función para graficarla
x <- seq(0, 2.5*Tam_Muestra, length = 10000)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
par(col.lab="#A32062", cex.lab=1.3, col.main="#FF4D00", cex.main=1.6,
     fg="#069900", mar=c(4,4,2,0.5), las=1)
hist(A, col = rainbow(8), density = 40, prob = T,
     main = bquote("HISTOGRAMA DEL COMPORTAMIENTO DE "~S^2),
     xlab = bquote(S^2), ylab = "Densidad")
lines(x, dchisq(x,df = Tam_Muestra-1), type = "l", col = "#D41243", lwd = 2)

ASINTÓTICAMENTE INSESGADOS

Recuerde que,

• \(S^2 = \frac{1}{n-1} \sum_{i=1}^{n}{(X_i - \overline{X})^2}\)
• \(S_n^2 = \frac{1}{n} \sum_{i=1}^{n}{(X_i - \overline{X})^2}\)

En el caso particular en que \(S_n^2\) denota la varianza de una muestra aleatoria de tamaño \(n\) de \(X \overset{d}{=} \mathcal{N}(\mu, \sigma^2)\) entonces, \[\frac{n S_n^2}{(n-1)} \overset{p}{\rightarrow} \sigma^2\]

Revisemos los estimadores \(S^2\) y \(S_n^2\) como estimadores de \(\sigma^2\). Aunque \(S_n^2\) resulta ser sesgado para \(\sigma^2\), el sesgo se hace pequeño cuando el tamaño muestral crece, pues,

\[\lim_{n\rightarrow\infty} B_{S_n^2} = \lim_{n\rightarrow\infty} -\frac{\sigma^2}{n} = 0\]

Para ilustrar los estimadores \(S^2\) y \(S_n^2\) en términos de estimación, podemos simular muestras provenientes de una distribución normal de media 5 y desviación estándar 5 con tamaños de muestra \(n = 2,10,25,\cdots ,1000\), y en cada muestra calculamos los dos estimadores. El comando en \(\texttt{R}\) y los resultados se muestran a continuación.

set.seed(1980)
n <- c(2,10,25,50,100,200,300,500,1000)
Var1 <- NULL
Var2 <- NULL

for (i in 1:length(n)) {
  Media = 5; Varianza = 25
  Data <- rnorm(n[i], mean = Media, sd = sqrt(Varianza))
  
  Var1[i] <- var(Data)
  Var2[i] <- (n[i]-1)/n[i]*var(Data)
}

par(col.lab="#00CC00", cex.lab=1.4, fg="#F67500", mar=c(4,4,0.5,0.5), las=1)
plot(Var1, type = "b", col = "#D81538", lwd = 2,
     ylim = c(min(Var2),max(Var1)), xaxt = "n",
     xlab = "Tamaño de muestra (# Valores simulados)",
     ylab = "Estimación de la varianza")
lines(Var2, type = "b", col = "#15D8B5", lwd = 2, pch = 2)
abline(h = Varianza, lwd = 2, col = "#FFA700")
axis(1, 1:length(n), n)

legend("bottomright", c("INSESGADO", "SESGADO"), lty = c(1,1), pch = c(1,2),
       col = c("#D81538","#15D8B5"), box.col = 'transparent', text.col = 1, bg = 'transparent')

Ilustración de las estimaciones de \(S^2\) y \(S_n^2\) como estimadores de \(\sigma^2\) en muestras provenientes de una \(\mathcal{N}(5,5^2)\).

Cómo podemos observar las estimaciones del estimador sesgado \(S_n^2\) siempre son inferiores que los del estimador insesgado \(S^2\) ; en segundo lugar, la diferencia entre los dos estimadores se hace cada vez más pequeña y los valores de ambos estimadores se acercan al parámetro teórico a medida que el tamaño muestral crece.

Por otro lado, aunque \(S_n^2\) subestima la varianza teórica, en la gráfica podemos observar momentos en que estuvo por encima de la varianza teórica, esto no es ninguna contradicción con el hecho de que \(S_n^2\) subestima a \(\sigma^2\), ya que el concepto de subestimación de un estimador indica que promediando todos los valores del estimador, da un valor inferior al parámetro, mas no indica que todas las estimaciones obtenidas son inferiores al parámetro.

PROBABILITY INTEGRAL TRANSFORM

Sea \(X\) una variable aleatoria con distribución \(F_X(x;\theta)\), la variable aleatoria \(U = F_X(x)\) tiene distribución \(U(0,1)\). De igual modo, si \(U\overset{d}{=}U(0,1)\) entonces \(X = F_X^{-1}(U)\) tiene distribución \(F_X(x;\theta)\). En resumen, tenemos que,

• \(X\overset{d}{=}F_X(x;\theta) {\color{Orange} \rightarrow } U=F_X(x;\theta)\overset{d}{=}U(0,1)\)
• \(U\overset{d}{=}U(0,1) {\color{DarkOrange} \rightarrow} X=F_X^{-1}(U)\overset{d}{=}F_X(x;\theta)\)

A continuación, comprobaremos dicho resultado mediante simulación, para ello usaremos la distribución normal, puede ser cualquier otra distribución, pero por costumbre usaremos esta. Primero generamos valores aleatorios de \(X \overset{d}{=} \mathcal{N}(\mu, \sigma^2)\) con la función rnorm(), seguidamente definiremos \(U = F_X(x)\) para lo cual tomaremos los valores de \(X\) y mediante pnorm() calcularemos \(F_X(x)\), finalmente graficaremos ambas funciones para observar lo propuesto por el teorema.

# Distribución de X
set.seed(4321)
par(col.axis="#EE2C2C", col.lab="#2B1A3E", col.main="#9E217C", mfrow=c(1,2), las=1)

mu <- 0; sigma <- 1
Min <- mu-3*sigma; Max <- mu+3*sigma
x1  <- seq(Min, Max, length = 1000)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
X <- rnorm(n = 10000, mean = mu, sd = sigma)
hist(X, prob = T, col = FullColor[5:15], xlim = c(Min,Max),
     main = "DISTRIBUCIÓN DE X", xlab = bquote(X%~%F[X](x,theta)), ylab = "")
lines(x1, dnorm(x1,mu,sigma), type = "l", col = "#CF0060", lwd = 2)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# Transformación U=F(X)
x1 <- seq(0, 1, length = 1000)
U <- pnorm(q = X, mean = mu, sd = sigma)
hist(U, prob = T, col = FullColor[5:15],
     main = bquote("DISTRIBUCIÓN DE "~U==F[X](x,theta)),
     xlab = bquote(U%~%U(0,1)), ylab = "")
lines(x1, dunif(x1,0,1), type = "l", col = "#CF0060", lwd = 2)

Para concluir realizaremos el procedimiento contrario, partiremos de simular valores de \(U\overset{d}{=}U(0,1)\) con la función runif() y realizaremos la transformación \(F_X^{-1}(U)\) mediante la función qnorm(), finalmente graficaremos ambas funciones para observar lo propuesto por el teorema.

# Distribución de U~U(0,1)
set.seed(1234)
par(mfrow=c(1,2))
par(col.axis="#37DA29", col.lab="#D40078", col.main="#FFB100", las=1)

U <- runif(n = 10000, min = 0, max = 1)
hist(U, prob=T, col = FullColorLight[c(1:5,15:19)],
     main = "DISTRIBUCIÓN DE U", xlab = expression(U%~%U(0,1)), ylab = "")
x1  <- seq(0, 1, length = 1000)
lines(x1, dunif(x1,0,1), type = "l", col = "#8CD52F", lwd = 2)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# Retransformación X=F^-1(U)
x1 <- seq(-3, 3, length = 1000)
X  <- qnorm(U, 0, 1)
hist(X, prob = T, col = FullColorLight[c(1:5,15:19)], xlim = c(-3,3),
     main = bquote("DISTRIBUCIÓN DE "~X=={F[X]^{-1}}(U)),
     xlab = bquote(X%~%F[X](x,theta)), ylab = "")
lines(x1, dnorm(x1,mu,sigma), type = "l", col = "#8CD52F", lwd = 2)

CANTIDAD PIVOTAL

INTERVALOS DE CONFIANZA PARA UNA MUESTRA

IC.Media_VarKnown <- function(Mean_Muestral, Varianza, n, NivConfianza) {
  Alpha = 1-NivConfianza
  Cuantil <- qnorm(1-(Alpha/2))
  Lim_Inf <- round(Mean_Muestral - (Cuantil*(sqrt(Varianza)/sqrt(n))),3)
  Lim_Sup <- round(Mean_Muestral + (Cuantil*(sqrt(Varianza)/sqrt(n))),3)
  cat(paste("I.C. DEL",NivConfianza*100,"% DE CONFIANZA =","\n    ",
            "(",Lim_Inf,",",Lim_Sup,")\n"))
}

IC.Media_VarUnknown <- function(Mean_Muestral, Var_Muestral, n, NivConfianza) {
  Alpha = 1-NivConfianza
  if(n<100) {
    Cuantil <- qt(1-(Alpha/2),n-1)
    Lim_Inf <- round(Mean_Muestral - (Cuantil*(sqrt(Var_Muestral)/sqrt(n))),3)
    Lim_Sup <- round(Mean_Muestral + (Cuantil*(sqrt(Var_Muestral)/sqrt(n))),3)
    cat(paste("I.C. DEL",NivConfianza*100,"% DE CONFIANZA =","\n    ",
              "(",Lim_Inf,",",Lim_Sup,")\n"))
  }
  else {
    Cuantil <- qnorm(1-(Alpha/2))
    Lim_Inf <- round(Mean_Muestral - (Cuantil*(sqrt(Var_Muestral)/sqrt(n))),3)
    Lim_Sup <- round(Mean_Muestral + (Cuantil*(sqrt(Var_Muestral)/sqrt(n))),3)
    cat(paste("I.C. DEL",NivConfianza*100,"% DE CONFIANZA =","\n    ",
              "(",Lim_Inf,",",Lim_Sup,")\n"))
  }
}

IC.Media_VarKnown(Mean_Muestral = 2.6, Varianza = 0.3^2, n = 36, NivConfianza = 0.95)

## I.C. DEL 95 % DE CONFIANZA = 
##      ( 2.502 , 2.698 )

IC.Media_VarKnown(Mean_Muestral = 2.6, Varianza = 0.3^2, n = 36, NivConfianza = 0.99)

## I.C. DEL 99 % DE CONFIANZA = 
##      ( 2.471 , 2.729 )

IC.Media_VarKnown(Mean_Muestral = 501, Varianza = 112^2, n = 500, NivConfianza = 0.99)

## I.C. DEL 99 % DE CONFIANZA = 
##      ( 488.098 , 513.902 )

Datos  <- c(2,3,3,5,1,1,2,1,1.5,2)
Media  <- mean(Datos)
Var    <- var(Datos)
Tamano <- length(Datos)
IC.Media_VarUnknown(Mean_Muestral = Media, Var_Muestral = Var, n = Tamano,
                    NivConfianza = 0.95)

## I.C. DEL 95 % DE CONFIANZA = 
##      ( 1.257 , 3.043 )

Datos  <- c(9.8,10.2,10.4,9.8,10,10.2,9.6)
Media  <- mean(Datos)
Var    <- var(Datos)
Tamano <- length(Datos)
IC.Media_VarUnknown(Mean_Muestral = Media, Var_Muestral = Var, n = Tamano,
                    NivConfianza = 0.95)

## I.C. DEL 95 % DE CONFIANZA = 
##      ( 9.738 , 10.262 )

IC.Var_Poblacional <- function(Var_Muestral, n, NivConfianza) {
  Alpha = 1-NivConfianza
  Cuantil1 <- qchisq(1-(Alpha/2),n-1)
  Cuantil2 <- qchisq(Alpha/2,n-1)
  Lim_Inf  <- round(((n-1)*Var_Muestral)/(Cuantil1),3)
  Lim_Sup  <- round(((n-1)*Var_Muestral)/(Cuantil2),3)
  cat(paste("I.C. DEL",NivConfianza*100,"% DE CONFIANZA =","\n    ",
            "(",Lim_Inf,",",Lim_Sup,")\n"))
}

IC.Var_Poblacional(Var_Muestral = 0.286, n = 10, NivConfianza = 0.95)

## I.C. DEL 95 % DE CONFIANZA = 
##      ( 0.135 , 0.953 )

IC.Proporcion <- function(Prop_Estimada, n, NivConfianza) {
  Alpha = 1-NivConfianza
  p = Prop_Estimada
  q = 1-Prop_Estimada
  Cuantil <- qnorm(1-(Alpha/2))
  Lim_Inf <- round(p - (Cuantil*sqrt((p*q/n))),3)
  Lim_Sup <- round(p + (Cuantil*sqrt((p*q/n))),3)
  cat(paste("I.C. DEL",NivConfianza*100,"% DE CONFIANZA =","\n    ",
            "(",Lim_Inf,",",Lim_Sup,")\n"))
}

IC.Proporcion(Prop_Estimada = 150/900, n = 900, NivConfianza = 0.95)

## I.C. DEL 95 % DE CONFIANZA = 
##      ( 0.142 , 0.191 )

prop.test(150, 900, conf.level = 0.95)

## 
##  1-sample proportions test with continuity correction
## 
## data:  150 out of 900, null probability 0.5
## X-squared = 398.67, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.1432251 0.1930061
## sample estimates:
##         p 
## 0.1666667

INTERVALOS DE CONFIANZA PARA DOS MUESTRAS

IC.Dif_Medias_VarsKnown <- function(Mean_1, Mean_2, Varianza1, Varianza2, n1, n2, NivConfianza) {
  Alpha = 1-NivConfianza
  Cuantil <- qnorm(1-(Alpha/2))
  Raiz    <- sqrt((Varianza1/n1)+(Varianza2/n2))
  Lim_Inf <- round((Mean_1-Mean_2) - (Cuantil*Raiz),3)
  Lim_Sup <- round((Mean_1-Mean_2) + (Cuantil*Raiz),3)
  cat(paste("I.C. DEL",NivConfianza*100,"% DE CONFIANZA =","\n    ",
            "(",Lim_Inf,",",Lim_Sup,")\n"))
}

IC.Dif_Medias_VarsIguales <- function(Mean_1, Mean_2, Var_Muestral1, Var_Muestral2, n1, n2, NivConfianza) {
  Alpha = 1-NivConfianza
  Cuantil <- qt(1-(Alpha/2),n1+n2-2)
  Raiz    <- sqrt((1/n1)+(1/n2))
  Sp      <- sqrt((((n1-1)*Var_Muestral1)+((n2-1)*Var_Muestral2))/(n1+n2-2))
  Lim_Inf <- round((Mean_1-Mean_2) - (Cuantil*Sp*Raiz),3)
  Lim_Sup <- round((Mean_1-Mean_2) + (Cuantil*Sp*Raiz),3)
  cat(paste("I.C. DEL",NivConfianza*100,"% DE CONFIANZA =","\n    ",
            "(",Lim_Inf,",",Lim_Sup,")\n"))
}

n1 <- 12; xbar <- 3.11; sx <- 0.771
n2 <- 10; ybar <- 2.04; sy <- 0.448

IC.Dif_Medias_VarsIguales(xbar,ybar,sx^2,sy^2,n1,n2,0.90)

## I.C. DEL 90 % DE CONFIANZA = 
##      ( 0.593 , 1.547 )

IC.Raz_Varianzas_MeanUnknown <- function(Var_Muestral1, Var_Muestral2, n1, n2, NivConfianza) {
  Alpha = 1-NivConfianza
  Cuantil_Inf <- qf(1-(Alpha/2), n1-1, n2-1)
  Cuantil_Sup <- qf(Alpha/2, n1-1, n2-1)
  Lim_Inf <- round(Var_Muestral1/(Var_Muestral2*Cuantil_Inf), 3)
  Lim_Sup <- round(Var_Muestral1/(Var_Muestral2*Cuantil_Sup), 3)
  cat(paste("I.C. DEL",NivConfianza*100,"% DE CONFIANZA =","\n    ",
            "(",Lim_Inf,",",Lim_Sup,")\n"))
}

n1 = 15; s1 = 3.07
n2 = 12; s2 = 0.8
IC.Raz_Varianzas_MeanUnknown(s1^2, s2^2, n1, n2, 0.98)

## I.C. DEL 98 % DE CONFIANZA = 
##      ( 3.43 , 56.903 )

PRUEBAS DE HIPÓTESIS

A continuación, realizaremos algunas pruebas de hipótesis con las funciones que tiene incorporadas \(\texttt{R}\), para lo cual estaremos bajo normalidad.

Prueba de hipótesis para una media (usando la distribución t de Student)

\[ \left\{ \begin{array}{ll} H_{0}: & \mu=10.5\\ H_{1}: & \mu<ó>ó\neq 10.5 \end{array} \right. \]

set.seed(2019)
mu    <- 10
sigma <- 5
x     <- rnorm(50, mu, sigma)

( A <- t.test(x, alternative = "two.sided", mu = 10.5, conf.level = 0.95) )

## 
##  One Sample t-test
## 
## data:  x
## t = -1.1968, df = 49, p-value = 0.2371
## alternative hypothesis: true mean is not equal to 10.5
## 95 percent confidence interval:
##   8.151332 11.095330
## sample estimates:
## mean of x 
##  9.623331

plot(A)

( B <- t.test(x, alternative = "less", mu = 10.5, conf.level = 0.95) )

## 
##  One Sample t-test
## 
## data:  x
## t = -1.1968, df = 49, p-value = 0.1186
## alternative hypothesis: true mean is less than 10.5
## 95 percent confidence interval:
##      -Inf 10.85139
## sample estimates:
## mean of x 
##  9.623331

plot(B)

( C <- t.test(x, alternative = "greater", mu = 10.5, conf.level = 0.95) )

## 
##  One Sample t-test
## 
## data:  x
## t = -1.1968, df = 49, p-value = 0.8814
## alternative hypothesis: true mean is greater than 10.5
## 95 percent confidence interval:
##  8.39527     Inf
## sample estimates:
## mean of x 
##  9.623331

plot(C)

Prueba de hipótesis para dos medias (usando la distribución t de Student)

set.seed(2019)
mu1     <- 10; mu2    <- 11
sigma1  <- 5; sigma2  <- 4
x     <- rnorm(50, mu1, sigma1)
y     <- rnorm(50, mu2, sigma2)

( A <- t.test(x, y, paired = T, alternative = "two.sided",
              mu= 1, var.equal = T, conf.level = 0.95,) )

## 
##  Paired t-test
## 
## data:  x and y
## t = -2.5493, df = 49, p-value = 0.01397
## alternative hypothesis: true difference in means is not equal to 1
## 95 percent confidence interval:
##  -2.7398954  0.5572322
## sample estimates:
## mean of the differences 
##               -1.091332

plot(A)

Prueba de hipótesis para dos varianzas (usando la distribución F de Fisher-Snedecor)

set.seed(2019)
mu    <- 10
sigma <- 5
x     <- rnorm(50, mu, sigma)

library(webr)
set.seed(2019)
mu1     <- 5; mu2    <- 10
sigma1  <- 1; sigma2 <- 2
x     <- rnorm(50, mu1, sigma1)
y     <- rnorm(50, mu2, sigma2)

( A <- var.test(x, y, paired = T, alternative = "two.sided", 
                ratio = 1, conf.level = 0.95) )

## 
##  F test to compare two variances
## 
## data:  x and y
## F = 0.4599, num df = 49, denom df = 49, p-value = 0.007539
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.2609840 0.8104357
## sample estimates:
## ratio of variances 
##           0.459903

plot(A)

POTENCIA \(\pi(\mu)\)

Para un sistema de hipótesis, puede haber más de una regla de decisión. En este caso, necesitamos comparar estas reglas de decisión en términos de la probabilidad de cometer los dos tipos de errores, para lo cual se define la función de potencia como,

\[P(\textit{Rechazar }H_0|H_0 \textit{ falsa})\]

Cabe resaltar que la función de potencia cambia su forma de calcularla y su gráfica dependiendo el sistema de hipótesis, si es sobre una o dos muestras y sobre que parámetro estamos realizando la prueba de hipótesis.

Ahora observemos un ejemplo de cómo calcular la función de potencia para la prueba de hipótesis para la media teórica, proveniente de una muestra bajo normalidad, y para el sistema de hipótesis,

\[H_0: \mu=\mu_0 \,\,\textit{ vs. }\,\, H_1: \mu\neq\mu_0\]

Tenemos que,

\[\begin{align*} \pi(\mu) &= P(\textit{Rechazar }H_0|H_0 \textit{ falsa}) \\ &= P(\bar{X}>\mu_0+\mathrm{Z}_{1-\frac{\alpha}{2}}\frac{\sigma}{\sqrt{n}})+P(\bar{X}<\mu_0-\mathrm{Z}_{1-\frac{\alpha}{2}}\frac{\sigma}{\sqrt{n}}) \\ &= P\left ( \frac{\bar{X}-\mu}{\sigma/\sqrt{n}}>\frac{\mu_0-\mu}{\sigma/\sqrt{n}}+\mathrm{Z}_{1-\frac{\alpha}{2}} \right )+P\left ( \frac{\bar{X}-\mu}{\sigma/\sqrt{n}}<\frac{\mu_0-\mu}{\sigma/\sqrt{n}}-\mathrm{Z}_{1-\frac{\alpha}{2}} \right ) \\ &= 1-\Phi\left ( \frac{\mu_0-\mu}{\sigma/\sqrt{n}}+\mathrm{Z}_{1-\frac{\alpha}{2}} \right )+\Phi\left ( \frac{\mu_0-\mu}{\sigma/\sqrt{n}}-\mathrm{Z}_{1-\frac{\alpha}{2}} \right ) \end{align*}\]

donde \(\Phi\) denota la función de distribución de una variable normal estándar.

Es claro que la función de potencia depende del nivel de significancia \(\alpha\) que se interpreta como la probabilidad de cometer el error Tipo I, y se determina de antemano. Como se sabe no se puede minimizar simultáneamente las probabilidades de cometer el error Tipo I y Tipo II (\(\alpha\) y \(\beta\)), entonces se espera que la función de potencia tome valores pequeños cuando \(\alpha\) es pequeño.

Los anteriores comentarios se pueden ilustrar graficando la función de potencia para diferentes tamaños de muestra. Utilizaremos la siguiente sintaxis y la gráfica se presenta a continuación.

Po_Norm <- function(mu, sigma, n, mu0, alpha) {
  1-pnorm((mu0-mu)*sqrt(n)/sigma+qnorm(1-alpha/2)) +
    pnorm((mu0-mu)*sqrt(n)/sigma-qnorm(1-alpha/2))
}

mu0   <- 0
sigma <- 1
alpha <- 0.05

n <- c(10, 30, 50)
Colores <- c("#1EBBD7", "#FC913A", "#FF4E50")

par(col.lab="#A32062", cex.lab=1.5, col.main="#FFB81C", cex.main=1.8,
    fg="#30A736", mar=c(5,5,3,0.5), las=1)
for(i in 1:length(n)){
  curve(Po_Norm(x,sigma,n[i],mu0,alpha), xlim = c(-2,2),
        add = i>1, xlab = bquote(mu), ylab = bquote(pi(mu)),
        bty = "n", col = Colores[i], lwd = 2.5, lty = i)
  legend(1.1, 0.05*i+0.1, legend = bquote(n==.(n[i])),
         col = Colores[i], lty = i, box.col = "transparent", bg = "transparent",
         text.col = "black", cex = 1.5)
}
title("COMPORTAMIENTO DE LA FUNCIÓN DE POTENCIA\n VARIANDO 'n'")

Por otro lado, también podemos graficar la función de potencia para diferentes niveles de significancia, para corroborar que entre mayor sea \(\alpha\), mayor es la potencia de la prueba. Una pequeña modificación al anterior código nos arroja la siguiente gráfica.

mu0   <- 0
sigma <- 1
n     <- 20

Alphas  <- c(0.03, 0.05, 0.1)
Colores <- c("#1EBBD7", "#FC913A", "#FF4E50")

par(col.lab="#A32062", cex.lab=1.5, col.main="#FFB81C", cex.main=1.7,
    fg="#30A736", mar=c(5,5,3,0.5), las=1)
for(i in 1:length(Alphas)){
  curve(Po_Norm(x,sigma,n,mu0,Alphas[i]), xlim = c(-2,2),
        add = i>1, xlab = bquote(mu), ylab = bquote(pi(mu)),
        bty = "n", col = Colores[i], lwd = 2.5, lty = i)
  legend(1.1, 0.05*i+0.1, legend = bquote(alpha==.(Alphas[i])),
         col = Colores[i], lty = i, box.col = "transparent", bg = "transparent",
         text.col = "black", cex = 1.5)
}
title(bquote("COMPORTAMIENTO DE LA FUNCIÓN DE POTENCIA VARIANDO "~alpha))

A continuación se expone el código \(\texttt{R}\) en el cual mediante simulación podemos apreciar que la potencia \(\pi(\mu)\), es decir, la \(P(\textit{Rechazar }H_0|H_0 \textit{ falsa})\) cuando \(\mu_{alterna}\) está cerca del verdadero de \(\mu\) es pequeña, pero aumenta rápidamente cuando \(\mu_{alterna}\) se aleja del valor verdadero.

mu    <- 150
sigma <- 10
n     <- 100
alpha <- 0.05

Valor_Critico  <- mu+qnorm(1-alpha/2,0,1)*(sigma/sqrt(n))
mu_Alterna     <- seq(from = 150, to = 155, by = 0.1)
Potencia <- NULL

for (i in 1:length(mu_Alterna)) {
  Cuantila    = (Valor_Critico-mu_Alterna[i])/(sigma/sqrt(n))
  Potencia[i] = 1-pnorm(Cuantila,0,1)
}

par(col.lab="#B6075D", cex.lab=1.5, col.main="#2EB82E", cex.main=1.8,
    fg="#7E00BF", mar=c(5,5,2,0.5), las=1)
plot(mu_Alterna, Potencia, type = "l", col = "#FF3300",
     main = "CURVA DE POTENCIA", lwd = 2,
     xlab = bquote(mu), ylab = bquote(pi(mu)))
abline(h = alpha, col = "#FFB81C", lty = 2)

A continuación realizaremos una pequeña modificación al grafico expuesto anteriormente representando está vez a \(\beta\) mediante la relación \(\beta(\mu)=1-\pi(\mu)\), así podemos apreciar que \(\beta(\mu)\), es decir, la \(P(\textit{No rechazar }H_0|H_0 \textit{ falsa})\) cuando \(\mu_{alterna}\) está cerca del verdadero de \(\mu\) es alta, pero disminuye rápidamente cuando \(\mu_{alterna}\) se aleja del valor verdadero.

par(col.lab="#001372", cex.lab=1.5, col.main="#03CEEB", cex.main=1.8,
    fg="#E32E01", mar=c(5,5,2,0.5), las=1)
plot(mu_Alterna, 1-Potencia, type = "l", col = "#37DA29",
     main = "CURVA DE BETA", lwd = 2,
     xlab = bquote(mu), ylab = bquote(beta(mu)))
abline(h = alpha, col = "#FFB81C", lty = 2)

TEST DE RAZÓN DE VEROSIMILITUD

Sea \(X_1, \cdots,X_n\) una muestra aleatoria de \(f_x(x;\theta)\). Sea

\[ \left\{ \begin{array}{ll} H_{0}: & \theta=\theta_0\\ H_{1}: & \theta=\theta_1 \end{array} \right. \]

un sistema de hipótesis simple, el test,

\[ \tau: \textit{Rechazar } H_0 \textit{ si } \lambda=\frac{L(\theta_0;x_1,\cdots,x_n)}{L(\theta_1;x_1,\cdots,x_n)}=\frac{\prod_{i=1}^{n}f_x(x_i;\theta_0)}{\prod_{i=1}^{n}f_x(x_i;\theta_1)}<\lambda_0 \]

recibe el nombre de test de razón de verosimilitud.

Ejemplo). Sea \(X \overset{d}{=} \textrm{P}(\lambda=3)\), se observa que \(x = 3\) y se quiere probar que,

\[ \left\{ \begin{array}{ll} H_{0}: & \theta=3.2\\ H_{1}: & \theta=4 \end{array} \right. \]

x   <- 3
Ho  <- 3.2
Ha  <- 4.0
( Lambda = dpois(x, Ho)/dpois(x, Ha) )

## [1] 1.139477

Observe que el numerador es mayor que el denominador, por lo tanto, podemos concluir que si la hipótesis nula es cierta \(\lambda\) tiende a ser “grande”.

Sea \(X \overset{d}{=} \textrm{P}(\lambda=4)\), se observa que \(x = 4\) y se quiere probar que,

\[ \left\{ \begin{array}{ll} H_{0}: & \theta=3\\ H_{1}: & \theta=3.8 \end{array} \right. \]

x   <- 4
Ho  <- 3
Ha  <- 3.8
( Lambda = dpois(x, Ho)/dpois(x, Ha) )

## [1] 0.8645422

Observe que el denominador es mayor que el numerador, por lo tanto, podemos concluir que si la hipótesis nula es falsa \(\lambda\) tiende a ser “pequeña”.

A continuación, haremos lo mismo, pero para muchos valores de \(x\), asumiendo \(H_0\) tanto cierta como falsa y suponiendo un tamaño de muestra \(n\).

Asumiendo \(H_0\) cierta,

Ho  <- 3.2
Ha  <- 4.0
n   <- 50
x   <- rpois(n, lambda = Ho)
Numerador   <- 1
Denominador <- 1
for (i in 1:n) {
  Numerador   = Numerador*dpois(x[i], Ho)
  Denominador = Denominador*dpois(x[i], Ha)
}
( Lambda = Numerador/Denominador )

## [1] 437.983

Asumiendo \(H_0\) falsa,

Ho  <- 3.2
Ha  <- 4.0
n   <- 50
x   <- rpois(n, lambda = Ha)
Numerador   <- 1
Denominador <- 1
for (i in 1:n) {
  Numerador   = Numerador*dpois(x[i], Ho)
  Denominador = Denominador*dpois(x[i], Ha)
}
( Lambda = Numerador/Denominador )

## [1] 0.009767348

Conclusión: Si la hipótesis nula es cierta \(\lambda\) tiende a ser “grande”, si por lo contrario, la hipótesis nula es falsa \(\lambda\) tiende a ser “pequeña”.

A continuación, realizaremos una generalización del ejemplo anterior, simulando (Tam_Simulacion = 100) muestras aleatorias de tamaño (n = 50) cada una, bajo \(H_0\) (suponiendo \(H_0\) verdadera) y bajo \(H_a\) (suponiendo \(H_0\) falsa).

set.seed(2020)
Ho  <- 3.2
Ha  <- 4.0
n   <- 50
Tam_Simulacion <- 100
Lambda_o <- NULL
for (i in 1:Tam_Simulacion) {
  x <- rpois(n, lambda = Ho)
  Numerador   <- 1
  Denominador <- 1
  for (j in 1:n) {
    Numerador   = Numerador*dpois(x[j], Ho)
    Denominador = Denominador*dpois(x[j], Ha)
  }
  Lambda_o[i]= Numerador/Denominador
}
summary(Lambda_o)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##     0.222    24.078   103.333  1235.542   727.120 30391.179

Intervalo_1   <- cut(x = Lambda_o, breaks = seq(0, 1000, 200), include.lowest = T)
Intervalo_Ho  <- cut(x = Lambda_o, breaks = seq(0, 100, 20), include.lowest = T)

par(col.axis="#AB0947", fg="#9FA09E", col.main="#AB0997", mfrow=c(1,2), las=2)
barplot(table(Intervalo_1), main = "Ho Verdadera", col = ColorG[2:6])
barplot(table(Intervalo_Ho), main = "Ho Verdadera", col = ColorG[7:11])

set.seed(2019)
Ho  <- 3.2
Ha  <- 4.0
n   <- 50
Tam_Simulacion <- 100
Lambda_a <- NULL
for (i in 1:Tam_Simulacion) {
  x   <- rpois(n, lambda = Ha)
  Numerador   <- 1
  Denominador <- 1
  for (j in 1:n) {
    Numerador   = Numerador*dpois(x[j], Ho)
    Denominador = Denominador*dpois(x[j], Ha)
  }
  Lambda_a[i]= Numerador/Denominador
}
summary(Lambda_a)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.000004 0.001946 0.023846 0.211265 0.096651 7.890001

Intervalo_1 <- cut(x = Lambda_a, breaks = seq(0, 100, 20), include.lowest=T)
Intervalo_Ha <- cut(x = Lambda_a, breaks = seq(0, 5, 0.5), include.lowest=T)

par(col.axis="#AB0947", fg="#9FA09E", col.main="#AB0997", mfrow=c(1,2), las=2)
barplot(table(Intervalo_1), main = "Ho Falsa", col = ColorG[2:6])
barplot(table(Intervalo_Ha), main = "Ho Falsa", col = ColorG[7:11])

Como observamos mediante el ejercicio de simulación, podemos percatarnos que efectivamente el comportamiento de \(\lambda\) es como los resultados teóricos proponen, a continuación, se expone la comparación del comportamiento de \(\lambda\).

par(mfrow=c(1,2), las=2)
barplot(table(Intervalo_Ho), col = rainbow(5,alpha = 0.7),
        main = bquote("Distribución de "~lambda~" bajo "~H[0]))
barplot(table(Intervalo_Ha), col = rainbow(5,alpha = 0.7),
        main = bquote("Distribución de "~lambda~" bajo "~H[a]))

Conclusión: Bajo \(H_0\) las lambdas son “grandes” y bajo \(H_a\) las lambdas son “pequeñas”. Es decir, debe rechazarse, de acuerdo con la simulación, la hipótesis nula si \(\lambda\) es “pequeña”.

TEST DE RAZÓN DE VEROSIMILITUD GENERALIZADA

Sea \(X_1, \cdots,X_nX_1,\cdots,X_n\) una muestra aleatoria de \(f_x(x;\theta)\). Sea

\[ \left\{ \begin{array}{ll} H_{0}: & \theta=\theta_0\\ H_{1}: & \theta<ó>ó\neq\theta_1 \end{array} \right. \]

un sistema de hipótesis compuesta, el test,

\[ \tau: \textit{Rechazar } H_0 \textit{ si } \lambda=\frac{\textit{Sup}_{\theta\in\Theta_0}L(\theta;x_1,\cdots,x_n)}{\textit{Sup}_{\theta\in\Theta}L(\theta;x_1,\cdots,x_n)}<\lambda_0 \]

recibe el nombre de test de razón de verosimilitud generalizada.

Recordemos que la función de verosimilitud es máxima en el EMV, para el siguiente ejemplo usaremos el caso de una Normal, pero puede ser cualquier otra distribución.

Suponga que se quiere probar,

\[ \left\{ \begin{array}{ll} H_{0}: & \mu=0 \\ H_{1}: & \mu>0 \end{array} \right. \]

A continuación, observaremos la distribución de \(\lambda\) bajo \(h_0\), mediante simulación, para lo cual tomaremos una muestra bajo \(h_0\) (\(\mu = 0\)), con lo cual deberíamos obtener que \(\lambda\) debería estar cercano a \(1\).

set.seed(2019)
n <- 10
Sigma <- 1
mu0 <- seq(from = -3, to = 0, length = 100)
mu1 <- seq(from = -3, to = 3, length = 100)
Lambda <- NULL
for (i in 1:100) {
  muHo <- 0
  x <- rnorm(n, muHo, Sigma)
  L_Num <- NULL
  L_Den <- NULL
  for(j in 1:100) {
    L_Num[j] = prod(dnorm(x,mu0[j],Sigma))
    L_Den[j] = prod(dnorm(x,mu1[j],Sigma))
  }
  Lambda[i] = max(L_Num)/max(L_Den)
}

( Lambda <- round(Lambda,2) )

##   [1] 1.00 1.00 0.90 1.00 0.90 1.00 1.00 0.42 1.00 1.00 1.00 1.00 1.00 0.16 0.95
##  [16] 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.52 1.00 1.00 1.00 1.00
##  [31] 1.00 1.00 1.00 1.00 1.00 1.00 0.99 1.00 0.48 0.87 1.00 1.00 0.40 1.00 1.00
##  [46] 1.00 0.99 0.27 0.95 1.00 1.00 0.99 0.95 1.00 0.76 0.86 1.00 0.31 1.00 1.00
##  [61] 1.00 1.00 1.00 1.00 1.00 1.00 0.15 1.00 1.00 1.00 1.00 0.52 1.00 1.00 0.37
##  [76] 1.00 0.52 1.00 0.55 1.00 0.85 0.23 1.00 1.00 1.00 1.00 0.59 1.00 0.98 1.00
##  [91] 1.00 1.00 0.91 1.00 0.96 0.12 0.73 0.80 1.00 1.00

par(col.axis="#CF0060", col.lab="#FF621E", col.main="#13A8FE", las=1)
hist(Lambda, col = rainbow(8), density = 40, xlim = c(0,1),
     main = bquote("Distribución Bajo "~H[0]), xlab = bquote(lambda), ylab = "Densidad")

Ahora observaremos la distribución de \(\lambda\) bajo \(h_a\), mediante simulación, para lo cual tomaremos una muestra bajo \(h_a\) (\(\mu>0\)), con lo cual deberíamos obtener que \(\lambda\) debería ser “pequeña”.

set.seed(2019)
n <- 10
Sigma <- 1
mu0 <- seq(from = -3, to = 0, length = 100)
mu1 <- seq(from = -3, to = 3, length = 100)
Lambda <- NULL
for (i in 1:100) {
  muHa <- runif(1, 0.1,3)
  x <- rnorm(n, muHa, Sigma)
  L_Num <- NULL
  L_Den <- NULL
  for(j in 1:100) {
    L_Num[j] = prod(dnorm(x,mu0[j],Sigma))
    L_Den[j] = prod(dnorm(x,mu1[j],Sigma))
  }
  Lambda[i] = max(L_Num)/max(L_Den)
}

( Lambda <- round(Lambda,2) )

##   [1] 0.00 0.00 0.00 0.02 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.01 0.09 0.01 0.00
##  [16] 1.00 0.00 0.00 1.00 0.92 0.00 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
##  [31] 0.01 0.15 0.00 0.02 0.09 0.00 0.40 0.00 1.00 0.00 0.00 0.01 0.00 0.19 0.00
##  [46] 0.00 0.00 0.02 0.00 0.09 0.14 0.00 0.02 0.01 0.00 1.00 0.55 0.05 0.00 0.32
##  [61] 0.10 0.00 0.00 0.00 0.00 0.52 0.00 0.00 0.00 0.00 0.91 0.00 1.00 0.93 0.01
##  [76] 0.00 0.01 0.00 0.00 0.01 0.00 0.00 0.00 0.38 0.00 0.00 0.00 0.03 0.00 0.27
##  [91] 0.00 0.00 0.05 0.00 0.00 0.00 0.01 0.00 0.36 0.00

par(col.axis="#CF0060", col.lab="#FF621E", col.main="#13A8FE", las=1)
hist(Lambda, col = rainbow(8), density = 40, xlim = c(0,1),
     main = bquote("Distribución Bajo "~H[a]), xlab = bquote(lambda), ylab = "Densidad")

RAZÓN DE VEROSIMILITUD MONÓTONA

Una familia de densidades \(f(x;\theta)\) tiene razón monótona de verosimilitudes en la estadística \(T=t(X_1, \cdots,X_n)\), si para dicha estadística,

\[\frac{L(\theta_1;x_1,\cdots,x_n)}{L(\theta_2;x_1,\cdots,x_n)}\]

Es función no creciente de \(t(x_1,\cdots,x_n)\) para \(\theta_1<\theta_2\), o es función no decreciente de \(t(x_1, \cdots,x_n)\) para \(\theta_1>\theta_2\).

Ejemplo 1). Distribución Poisson <\(X \overset{d}{=} \textrm{P}(\lambda)\)>.

\[\begin{align*} f(x;\theta) &= \frac{e^{-\theta}\theta^{x}}{x!} \\ L(\theta;x_1,\cdots,x_n) &= \frac{e^{-n\theta}\theta^{\sum_{i=1}^{n}X_i}}{\prod_{i=1}^{n}X_i!} \\ \frac{L(\theta_1;x_1,\cdots,x_n)}{L(\theta_2;x_1,\cdots,x_n)} &= e^{n(\theta_1-\theta_2)}\left ( \frac{\theta_1}{\theta_2} \right )^{\sum_{i=1}^{n}X_i} = y \end{align*}\]

Theta <- 2; Theta1 <- 3; Theta2 <- 4
n <- 20

Lambda <- NULL
Suma   <- NULL
for (i in 1:50) {
  x <- rpois(n, lambda = Theta)
  Suma[i] = sum(x)
  Lambda[i] = exp(n*(Theta2-Theta1))*(Theta1/Theta2)^(Suma[i]) 
}

par(col.axis="#DB165B", col.lab="#FF1414", cex.lab=1.5, col.main="#FF621E", cex.main=1.8,
    fg="#FF8D00", mar=c(5,5,2,0.5))
plot(Suma, Lambda, type = "p", pch = 18, col = 3, cex = 2,
     main = bquote(X[i]~"~ Poisson ("~lambda~") ES MONÓTONA DECRECIENTE"),
     xlab = bquote(sum(x[i])), ylab = bquote(lambda))

La anterior función es decreciente (no creciente) en \(\sum_{i=1}^{n}X_i\) puesto que \(\theta_1<\theta_2\), además, cuando \(\sum_{i=1}^{n}X_i\) aumenta \(y\) disminuye. Por lo tanto, la distribución Poisson tiene razón de verosimilitud monótona en \(\sum_{i=1}^{n}X_i\).

Ejemplo 2). Distribución Exponencial <\(X \overset{d}{=} \textrm{Exp}(\theta)\)>.

\[\begin{align*} f(x;\theta) &= \theta e^{-\theta x} \\ L(\theta;x_1,\cdots,x_n) &= \theta^{n} e^{-\theta \sum_{i=1}^{n}X_i} \\ \frac{L(\theta_1;x_1,\cdots,x_n)}{L(\theta_2;x_1,\cdots,x_n)} &= \left ( \frac{\theta_1}{\theta_2} \right )^{n} e^{-(\theta_1-\theta_2)\sum_{i=1}^{n}X_i} = y \end{align*}\]

Theta <- 2; Theta1 <- 3; Theta2 <- 4
n <- 20

Lambda <- NULL
Suma   <- NULL
for (i in 1:50) {
  x <- rexp(n, rate = Theta)
  Suma[i] = sum(x)
  Lambda[i] = (Theta1/Theta2)^n*exp((Theta2-Theta1)*Suma[i]) 
}

par(col.axis="#DB165B", col.lab="#FF1414", cex.lab=1.5, col.main="#FF621E", cex.main=1.8,
    fg="#FF8D00", mar=c(5,5,2,0.5))
plot(Suma, Lambda, type = "p", pch = 18, col = 4, cex = 2,
     main = bquote(X[i]~"~ Exp ("~theta~") ES MONÓTONA CRECIENTE"),
     xlab = bquote(sum(x[i])), ylab = bquote(theta))

La anterior función es creciente (no decreciente) en \(\sum_{i=1}^{n}X_i\). Luego la distribución Exponencial tiene razón de verosimilitud monótona en \(\sum_{i=1}^{n}X_i\).

Este obra cuyo autor es Jeison Mauricio Alarcón Becerra está bajo una licencia de Reconocimiento-NoComercial-SinObraDerivada 4.0 Internacional de Creative Commons.
Creado a partir de la obra en https://rpubs.com/JeisonAlarcon.