Para demostrar la versatilidad de esta técnica, generaré un conjunto de datos de 6 dimensiones y cada columna contendrá datos numéricos de varias distribuciones aleatorias. Aquí están las distribuciones y parámetros específicos de cada uno de los que usaré:
X1 ~ distribución t (v = 3)
X2 ~ Binomio (norte = 10, pag = 0,4)
X3 ~ Gama (α = 3, b = 2)
X4 ~ Uniforme (un = 0, b = 1)
X5 ~ Exponencial (λ = 50)
X6 ~ Poison (λ = 2)
Tenga en cuenta que cada columna de este conjunto de datos contiene variables independientes, ya que ninguna columna debe correlacionarse con otra ya que se crearon de forma independiente. Este no es un requisito para utilizar este método. Se hizo de esta manera simplemente por simplicidad y para demostrar la naturaleza paradójica de este resultado. Si no está completamente familiarizado con alguna o todas estas distribuciones, incluiré una imagen simple de cada una de las columnas univariadas de los datos generados aleatoriamente. Esta es simplemente una iteración de 1000 variables aleatorias generadas a partir de cada una de las distribuciones antes mencionadas.
De los histogramas anteriores debe quedar claro que no todas estas variables siguen una distribución normal, lo que implica que el conjunto de datos en su conjunto no es normal multivariado.
Como se conocen las distribuciones verdaderas de cada una, conocemos los promedios verdaderos de cada una. El promedio de este conjunto de datos multivariado se puede expresar en forma vectorial y cada entrada de fila representa el promedio de la variable respectivamente. En este ejemplo,
Conocer los promedios verdaderos de cada variable nos permitirá medir qué tan cerca está la media muestral, o el estimador de James Stein implica que cuanto más cerca, mejor. A continuación se muestra el experimento que ejecuté en código R que generó cada una de las 6 variables aleatorias y las probé con los promedios reales utilizando el error cuadrático medio. Luego, este experimento se realizó 10.000 veces utilizando cuatro tamaños de muestra diferentes: 5, 50, 500 y 5.000.
set.seed(42)
## Function to calculate Mean Squared Error ##
mse <- function(x, true_value)
return( mean( (x - true_value)^2 ) )
## True Average ##
mu <- c(0, 4, 1.5, 0.5, 0.02, 2)
## Store Average and J.S. Estimator Errors ##
Xbar.MSE <- list(); JS.MSE <- list()
for(n in c(5, 50, 500, 5000)){ # Testing sample sizes of 5, 30, 200, and 5,000
for(i in 1:1e4){ # Performing 10,000 iterations## Six Random Variables ##
X1 <- rt(n, df = 3)
X2 <- rbinom(n, size = 10, prob = 0.4)
X3 <- rgamma(n, shape = 3, rate = 2)
X4 <- runif(n)
X5 <- rexp(n, rate = 50)
X6 <- rpois(n, lambda = 2)
X <- cbind(X1, X2, X3, X4, X5, X6)
## Estimating Std. Dev. of Each and Standardizing Data ##
sigma <- apply(X, MARGIN = 2, FUN = sd)
## Sample Mean ##
Xbar <- colMeans(X)
## J.S. Estimator ##
JS.Xbar <- james_stein_estimator(Xbar=Xbar, sigma2=sigma/n)
Xbar.MSE[[as.character(n)]][i] <- mse(Xbar, mu)
JS.MSE[[as.character(n)]][i] <- mse(JS.Xbar, mu)
}
}
sapply(Xbar.MSE, mean) # Avg. Sample Mean MSE
sapply(JS.MSE, mean) # Avg. James-Stein MSE
De los 40.000 senderos, el MSE promedio total de cada tamaño de muestra se calcula ejecutando las dos últimas líneas. Los resultados de cada uno se pueden ver en la siguiente tabla.
Los resultados de esta simulación muestran que el estimador de James-Stein es consistentemente mejor que la media muestral usando el MSE, pero que esta diferencia disminuye a medida que aumenta el tamaño de la muestra.