#install.packages("kableExtra")
#install.packages("moments")
#install.packages("gganimate")
#install.packages("gifski")
#install.packages("ggplot2")
#install.packages("magick")
#install.packages("reshape2")
#install.packages("plotly")
La simulación ayuda a entender y validar las propiedades de los estimadores estadísticos como son: insesgadez, eficiencia y consistencia, principalmente. El siguiente problema permite evidenciar las principales características de un grupo de estimadores propuestos para la estimación de un parámetro asociado a un modelo de probabilidad.
Sean \(X_{1}\), \(X_{2}\), \(X_{3}\) y \(X_{4}\), una muestra aleatoria de tamaño \(n=4\) cuya población la conforma una distribución exponencial con parámetro \(\theta\) desconocido. Determine las características de cada uno de los siguientes estimadores propuestos:
\[\begin{aligned} \hat{\theta}_{1}=\frac{X_{1}+X_{2}}{6}+\frac{X_{3}+X_{4}}{3} \end{aligned}\]
\[\begin{aligned} \hat{\theta}_{2}=\frac{(X_{1}+2X_{2}+3X_{3}+4X_{4})}{5} \end{aligned}\]
\[\begin{aligned} \hat{\theta}_{3}=\frac{X_{1}+X_{2}+X_{3}+X_{4}}{4} \end{aligned}\]
\[\begin{aligned} \hat{\theta}_{4}=\frac{\min(X_{1}+X_{2}+X_{3}+X_{4})+\max(X_{1}+X_{2}+X_{3}+X_{4})}{2} \end{aligned}\]
Se consideran los siguientes conceptos:
Insesgadez: Un estimador \(\hat{\theta}\) es insesgado si su valor esperado es igual al valor del parámetro verdadero \(\theta\).
\[ \begin{aligned} E(\hat{\theta})=\theta \end{aligned} \]
\(\hat{\theta}_1 = \frac{X_1 + X_2}{6} + \frac{X_3 + X_4}{3}\)
\(E(\hat{\theta}_1) = E(\frac{X_1 + X_2}{6} + \frac{X_3 + X_4}{3}) = \frac {E(X_1 + X_2)}{6} + \frac{E(X_3 + X_4)}{3}\)
\(E(\hat{\theta}_1) = \frac {E(X_1) + E(X_2)}{6} + \frac{E(X_3) + E(X_4)}{3}\)
\(E(\hat{\theta}_1) = \frac {{\theta}+{\theta}}{6} + \frac{{\theta} + {\theta}}{3} = \frac {2{\theta}}{6} + \frac{2{\theta}}{3}\)
\(E(\hat{\theta}_1) = \frac {{\theta}}{3} + \frac{2{\theta}}{3} = \frac{3{\theta}}{3}\)
\(E(\hat{\theta}_1) = {\theta}\)
Por tanto, se comprueba que el estimador 1 es insesgado.
\(\hat{\theta}_2 = \frac{X_1 + 2X_2 + 3X_3 + 4X_4}{5}\)
\(E(\hat{\theta}_2) = E(\frac{X_1 + 2X_2 + 3X_3 + 4X_4}{5}) = \frac{E(X_1 + 2X_2 + 3X_3 + 4X_4)}{5}\)
\(E(\hat{\theta}_2) = \frac{E(X_1) + 2E(X_2) + 3E(X_3) + 4E(X_4)}{5}\)
\(E(\hat{\theta}_2) = \frac {{\theta}+2{\theta}+3{\theta}+4{\theta}}{5} = \frac {10{\theta}}{5}\)
\(E(\hat{\theta}_2) = 2{\theta}\)
Por tanto, se comprueba que el estimador 2 no es insesgado.
\(\hat{\theta}_3 = \frac{X_1 + X_2 + X_3 + X_4}{4}\)
\(E(\hat{\theta}_3) = E(\frac{X_1 + X_2 + X_3 + X_4}{4}) = \frac{E(X_1 + X_2 + X_3 + X_4)}{4}\)
\(E(\hat{\theta}_3) = \frac{E(X_1) + E(X_2) + E(X_3) + E(X_4)}{4}\)
\(E(\hat{\theta}_3) = \frac {{\theta}+{\theta}+{\theta}+{\theta}}{5} = \frac {4{\theta}}{4}\)
\(E(\hat{\theta}_2) = {\theta}\)
Por tanto, se comprueba que el estimador 3 es insesgado.
\(\hat{\theta}_4 = \frac{min\{X_1, X_2, X_3, X_4\}+max\{X_1, X_2, X_3, X_4\}}{2}\)
\(E(\hat{\theta}_4) = E(\frac{min\{X_1, X_2, X_3, X_4\}+max\{X_1, X_2, X_3, X_4\}}{2})\)
Debido a la complejidad analitica que plantea la función mínimo y máximo, la insesgadez se analizará a través de una simulación.
Eficiencia: Sea \(\hat{\theta}\) un estimador de un parámetro \({\theta}\), y sea \(Var(\hat{\theta})\) su varianza. El estimador \(\hat{\theta}\) es eficiente si su varianza es la menor entre todas las varianzas de los estimadores no sesgados de \({\theta}\). Es decir:
\[ \begin{aligned} Var(\hat{\theta}) \leq Var(\tilde{\theta}) \quad \forall \, \tilde{\theta} \, \text{no sesgado}. \end{aligned} \]
\(\hat{\theta}_1 = \frac{X_1 + X_2}{6} + \frac{X_3 + X_4}{3}\)
\(Var(\hat{\theta}_1) = Var(\frac{X_1 + X_2}{6} + \frac{X_3 + X_4}{3}) = \frac{Var(X_1 + X_2)}{6^2} + \frac{Var(X_3 + X_4)}{3^2}\)
\(Var(\hat{\theta}_1) = \frac{Var(X_1) + Var(X_2)}{36} + \frac{Var(X_3) + Var(X_4)}{9} = \frac{\sigma^2 + \sigma^2}{36} + \frac{\sigma^2 + \sigma^2}{9}\)
\(Var(\hat{\theta}_1) = \frac{2\sigma^2}{36} + \frac{2\sigma^2}{9} = \frac{\sigma^2}{18} + \frac{4\sigma^2}{18}\)
\(Var(\hat{\theta}_1) = \frac{5\sigma^2}{18}\)
\(\hat{\theta}_2 = \frac{X_1 + 2X_2 + 3X_3 + 4X_4}{5}\)
\(Var(\hat{\theta}_2) = Var(\frac{X_1 + 2X_2 + 3X_3 + 4X_4}{5}) = \frac{Var(X_1 + 2X_2 + 3X_3 + 4X_4)}{5}\)
\(Var(\hat{\theta}_2) = \frac{Var(X_1) + 2^2Var(X_2) + 3^2Var(X_3) + 4^2Var(X_4)}{5^2}\)
\(Var(\hat{\theta}_2) = \frac {\sigma^2+2\sigma^2+3\sigma^2+4\sigma^2}{25} = \frac {10\sigma^2}{25}\)
\(Var(\hat{\theta}_2) = \frac {2{\sigma^2}}{5}\)
\(\hat{\theta}_3 = \frac{X_1 + X_2 + X_3 + X_4}{4}\)
\(Var(\hat{\theta}_3) = Var(\frac{X_1 + X_2 + X_3 + X_4}{4}) = \frac{Var(X_1 + X_2 + X_3 + X_4)}{4}\)
\(Var(\hat{\theta}_3) = \frac{Var(X_1) + Var(X_2) + Var(X_3) + Var(X_4)}{4}\)
\(Var(\hat{\theta}_3) = \frac {\sigma^2+\sigma^2+\sigma^2+\sigma^2}{4} = \frac {4\sigma^2}{4}\)
\(Var(\hat{\theta}_2) = {\sigma^2}\)
\(\hat{\theta}_4 = \frac{min\{X_1, X_2, X_3, X_4\}+max\{X_1, X_2, X_3, X_4\}}{2}\)
\(Var(\hat{\theta}_4) = Var(\frac{min\{X_1, X_2, X_3, X_4\}+max\{X_1, X_2, X_3, X_4\}}{2})\)
Debido a la complejidad analitica que plantea la función mínimo y máximo, la eficiencia se analizará a través de una simulación.
Consistencia: Un estimador \(\hat{\theta}\) es consistente si, a medida que el tamaño de la muestra aumenta, el estimador converge al valor del parámetro verdadero \(\theta\).
Se generan las variables \(X_{1}\), \(X_{2}\), \(X_{3}\) y \(X_{4}\) y se calculan los estimadores \(\hat{\theta}_{1}\), \(\hat{\theta}_{2}\), \(\hat{\theta}_{3}\) y \(\hat{\theta}_{4}\), considerando la distribución exponencial con \(\theta\)\(=100\), para tamaños de muestra de \(n=20, 50, 100, 1000\). Luego, se evalua la eficiencia y la insesgadez para cada estimador y tamaño de muestra.
library(kableExtra)
set.seed(999)
n=20
theta=100
x1=rexp(n,1/theta)
x2=rexp(n,1/theta)
x3=rexp(n,1/theta)
x4=rexp(n,1/theta)
data.X <- data.frame(x1,x2,x3,x4)
theta1=(x1+x2)/6 + (x3+x4)/3
theta2=(x1 + 2*x2 + 3*x3 + 4*x4)/5
theta3=(x1 + x2 + x3 + x4)/4
min.T<-apply(data.X,1,min)
max.T<-apply(data.X,1,max)
theta4=(min.T + max.T)/2
data.T20 = data.frame(theta1,theta2,theta3,theta4)
resultados20 <- data.frame(n,
Media=round(apply(data.T20,2,mean),3),
Desviacion=round(apply(data.T20,2,sd),3),
Varianza=round(apply(data.T20,2,var),3),
Insesgadez=round(apply(data.T20,2,mean)-theta,3),
Consistencia=round((apply(data.T20,2,mean)-theta)/theta,3)
)
kable(resultados20, "html", caption = "Resultados de estimadores para n = 20 y θ = 100") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE)
| n | Media | Desviacion | Varianza | Insesgadez | Consistencia | |
|---|---|---|---|---|---|---|
| theta1 | 20 | 103.755 | 40.146 | 1611.697 | 3.755 | 0.038 |
| theta2 | 20 | 210.000 | 85.333 | 7281.644 | 110.000 | 1.100 |
| theta3 | 20 | 101.702 | 42.143 | 1776.011 | 1.702 | 0.017 |
| theta4 | 20 | 122.472 | 54.696 | 2991.612 | 22.472 | 0.225 |
set.seed(999)
n=50
theta=100
x1=rexp(n,1/theta)
x2=rexp(n,1/theta)
x3=rexp(n,1/theta)
x4=rexp(n,1/theta)
data.X <- data.frame(x1,x2,x3,x4)
theta1=(x1+x2)/6 + (x3+x4)/3
theta2=(x1 + 2*x2 + 3*x3 + 4*x4)/5
theta3=(x1 + x2 + x3 + x4)/4
min.T<-apply(data.X,1,min)
max.T<-apply(data.X,1,max)
theta4=(min.T + max.T)/2
data.T50 = data.frame(theta1,theta2,theta3,theta4)
resultados50 <- data.frame(n,
Media=round(apply(data.T50,2,mean),3),
Desviacion=round(apply(data.T50,2,sd),3),
Varianza=round(apply(data.T50,2,var),3),
Insesgadez=round(apply(data.T50,2,mean)-theta,3),
Consistencia=round((apply(data.T50,2,mean)-theta)/theta,3)
)
kable(resultados50, "html", caption = "Resultados de estimadores para n = 50 y θ = 100") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE)
| n | Media | Desviacion | Varianza | Insesgadez | Consistencia | |
|---|---|---|---|---|---|---|
| theta1 | 50 | 103.820 | 40.872 | 1670.510 | 3.820 | 0.038 |
| theta2 | 50 | 207.514 | 81.086 | 6574.978 | 107.514 | 1.075 |
| theta3 | 50 | 104.051 | 38.840 | 1508.550 | 4.051 | 0.041 |
| theta4 | 50 | 119.938 | 52.507 | 2756.968 | 19.938 | 0.199 |
set.seed(999)
n=100
theta=100
x1=rexp(n,1/theta)
x2=rexp(n,1/theta)
x3=rexp(n,1/theta)
x4=rexp(n,1/theta)
data.X <- data.frame(x1,x2,x3,x4)
theta1=(x1+x2)/6 + (x3+x4)/3
theta2=(x1 + 2*x2 + 3*x3 + 4*x4)/5
theta3=(x1 + x2 + x3 + x4)/4
min.T<-apply(data.X,1,min)
max.T<-apply(data.X,1,max)
theta4=(min.T + max.T)/2
data.T100 = data.frame(theta1,theta2,theta3,theta4)
resultados100 <- data.frame(n,
Media=round(apply(data.T100,2,mean),3),
Desviacion=round(apply(data.T100,2,sd),3),
Varianza=round(apply(data.T100,2,var),3),
Insesgadez=round(apply(data.T100,2,mean)-theta,3),
Consistencia=round((apply(data.T100,2,mean)-theta)/theta,3)
)
kable(resultados100, "html", caption = "Resultados de estimadores para n = 100 y θ = 100") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE)
| n | Media | Desviacion | Varianza | Insesgadez | Consistencia | |
|---|---|---|---|---|---|---|
| theta1 | 100 | 100.862 | 45.743 | 2092.461 | 0.862 | 0.009 |
| theta2 | 100 | 199.527 | 92.161 | 8493.582 | 99.527 | 0.995 |
| theta3 | 100 | 101.659 | 42.660 | 1819.852 | 1.659 | 0.017 |
| theta4 | 100 | 118.177 | 52.568 | 2763.421 | 18.177 | 0.182 |
set.seed(999)
n=1000
theta=100
x1=rexp(n,1/theta)
x2=rexp(n,1/theta)
x3=rexp(n,1/theta)
x4=rexp(n,1/theta)
data.X <- data.frame(x1,x2,x3,x4)
theta1=(x1+x2)/6 + (x3+x4)/3
theta2=(x1 + 2*x2 + 3*x3 + 4*x4)/5
theta3=(x1 + x2 + x3 + x4)/4
min.T<-apply(data.X,1,min)
max.T<-apply(data.X,1,max)
theta4=(min.T + max.T)/2
data.T1000 = data.frame(theta1,theta2,theta3,theta4)
resultados1000 <- data.frame(n,
Media=round(apply(data.T1000,2,mean),3),
Desviacion=round(apply(data.T1000,2,sd),3),
Varianza=round(apply(data.T1000,2,var),3),
Insesgadez=round(apply(data.T1000,2,mean)-theta,3),
Consistencia=round((apply(data.T1000,2,mean)-theta)/theta,3)
)
kable(resultados1000, "html", caption = "Resultados de estimadores para n = 1.000 y θ = 100") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE)
| n | Media | Desviacion | Varianza | Insesgadez | Consistencia | |
|---|---|---|---|---|---|---|
| theta1 | 1000 | 98.680 | 51.575 | 2659.979 | -1.320 | -0.013 |
| theta2 | 1000 | 196.302 | 105.304 | 11088.885 | 96.302 | 0.963 |
| theta3 | 1000 | 98.836 | 49.034 | 2404.321 | -1.164 | -0.012 |
| theta4 | 1000 | 116.327 | 61.567 | 3790.462 | 16.327 | 0.163 |
par(mfrow = c(2, 2))
boxplot(data.T20,main = "n=20")
abline(h=theta, col="red", las=1)
grid()
boxplot(data.T50,main = "n=50")
abline(h=theta, col="red", las=1)
grid()
boxplot(data.T100,main = "n=100")
abline(h=theta, col="red", las=1)
grid()
boxplot(data.T1000,main = "n=1.000")
abline(h=theta, col="red", las=1)
grid()
En relacion con la insesgadez, los parametros \(\hat{\theta}_{1}\) y \(\hat{\theta}_{3}\) demuestran ser insesgados ya que el valor esperado de dichos estimadores es igual a \({\theta}\). Este hecho también se comprueba en la simulación al definir un \({\theta} = 100\) y evidenciar que dichos estimadores son los más cercanos a \(100\). Además, logran demostrar consistencia al converger al parámetro verdadero de \({\theta}\) a medida que aumenta el tamaño de la muestra.
Con respecto a la eficiencia, se considera que \(\hat{\theta}_{3}\) es el estimador más eficiente, ya que en la demostración \(Var(\hat{\theta}_3) = {\theta}\) y en la simulación dicho estimador presenta la menor varianza.
Para el resto de estimadores, \(\hat{\theta}_{2}\) es el de menor insesgadez al encontrarse muy lejos del parámetro \(\theta=100\), duplicando su valor; demuestra menor eficiencia dado que tiene la mayor varianza entre los cuatro estimadores; y por último es el de menor consistencia ya que, aunque el tamaño de la muestra aumente, no converge a \(\theta\).
El estimador \(\hat{\theta}_{4}\) tampoco cumple los criterios de manera óptima, al no lograr ser insesgado, no presentar la menor varianza y no ser lo suficiente consistente, al menos para los tamaños de muestra considerados.
En conclusión, el estimador \(\hat{\theta}_{3}\) es el más apropiado.