Simulación del efecto Stein

Mostrando el efecto Shrinkage

Caso normal

El estimador de James-Stein es un tipo de estimador que contrae el vector de medias hacia el origen tomado, con mayor contracción conforme la norma de $X\sim\mathcal{N}_p(\theta, \sigma^2 I_p)$ disminuye. $$\hat \delta_{𝐽𝑆}=\left(1−\frac{(𝑝−2)𝜎^2}{‖X‖^2}\right)X. $$

Iniciamos las simulaciones suponiendo una normal Gaussianana p-multivariante, fijando $p=20$, para ver la contracción que se produce en cada componente hacia el origen tomado.

Primero vemos el caso en el que el término multiplicativo del estimador de James-Stein es completamente positivo.

Otra casuística se da cuando el término multiplicativo es negativo y por lo tanto, obtenemos un mal resultado:

A continuación, vemos el caso en que no todas las componentes del estimador JS son menores a las del EMV, pero que aún así hay una contracción significativa.

Mostrando el ECM en función de la dimensión

Simulamos $N=50000$ muestras de un vector aleatorio p-variante normal $X\sim \mathcal{N}_p(\theta, \sigma^2 I_p)$, con la dimensión $p\geq 3$ variable. Es de esperar que el valor observado $X$ tenga un ECM igual a $p\sigma^2$.

En el caso Gaussiano, vemos que el estimador de James-Stein cancela completamente el aumento lineal del ECM del valor observado (estimador usual) con un valor constante de $2,5$ aproximadamente, e incluso algo mejor al tomar la versión positiva del estimador.

Si el encogimiento lo hacemos hacia el vector de medias $\theta=(1, ..., 1)$ seguimos teniendo cierta mejora aunque no tan significativa.

Vemos también que al tomar $n\geq1$ observaciones de cada muestra, el estimador de James-Stein (JS) se obtiene a partir de $$\hat \delta_{𝐽𝑆}=\left(1−\frac{(𝑝−2)𝜎^2}{n‖\bar X‖^2}\right)\bar X. $$ Donde $\bar X$ es la matriz de medias de dimensión $N\times p$ de cada $n$.

Simulamos el caso $n=10$. El ECM esperado de $X$ será pues $\dfrac{p\sigma^2}{n}=\dfrac{p\sigma^2}{10}$. De este modo, el caso Gaussiano anterior queda graficado de la siguiente forma:

Es decir, los estimadores JS y JS+ siguen cancelando el aumento del ECM de $X$ pero esta vez con un valor constante muy cercano a $0$. Esto muestra que cuanta más información no dispar tenemos, menor es el riesgo.

Cabe remarcar que, el estimador JS centrado (y derivados) no deja de ser el estimador empírico Bayesiano de James-Stein, es cual calculamos de la siguiente manera: $$\hat \delta_{𝐽𝑆}=\bar X + \left(1−\frac{(𝑝−3)\sigma^2}{n‖X -\bar X‖^2}\right)(X - \bar X). $$

Tomando ahora $p\geq4$.

Suponemos entonces, que nuestro vector de medias sigue una distribución previa, $\theta\sim (\mu, \tau^2)$ con $\mu=0$ y $\tau$=1. Y vemos el caso $n=10$.

Y si $n=1$, comparándolo con el caso $X\sim\mathcal(\theta,I)$ con $\theta$ Gaussiana, vemos que al añadir varianza al vector de medias, los estimadores presentan más desviación.

Pero si $\tau=5$ (aumentamos la varianza) entonces reducimos la mejora:

Es decir, cuanta mayor sea la varianza del vector de medias más se aproximan nuestros estimadores al EMV, y por lo tanto, menos mejora. Además, comparando los casos $n=10$ y $n=1$ de $X\sim \mathcal{N}(\theta, I_p)$, $\theta\sim\mathcal{N}(1, 1)$ observamos que al tener más observaciones con varianza de $\theta$ obtenemos menor mejora.

Mostrando el riesgo en función de la norma del vector de medias

Creamos una función para graficar los ECM según los distintos tipos de estimadores: el estimador usual (valor observado), el de James-Stein, el centrado y el de la parte positiva centrado para ver las mejoras que aportan según el tipo de vector de medias que tenemos. Fijamos $Var(X)=\sigma^2=1$.

Dicha función la calculamos en base al vector de medias así como en base a las observaciones que tengamos de cada muestra, $n\geq1$.

Veremos que cuanto vayor varianza del vector de medias, menor mejora tendremos (Efron, 1992). Es decir, cuando las coordenadas de $\theta$ tengan magnitud similar será cuando obtendremos mejores ganancias (menor riesgo).

Simulamos primero cuando la media tiene una distribución previa, $\theta\sim \mathcal{N}(\mu, \tau^2)$. Observamos el caso $\mu=1$ y $\tau=1$.

Si aumentamos la varianza, $\theta\sim \mathcal{N}(1, 5^2)$

Con $\theta$ Gaussiana:

Ahora si $\theta\sim Pois(\lambda)$, $\lambda=1$.

Cuando la media es cercana a un múltiplo de $\theta=(1, ..., 1)$

Cuando la media es cercana a un múltiplo de $\theta=(1, 2, ..., p)$

Cuando la media es cercana a un múltiplo de $\theta=(-1, ...,-1, 1, ..., 1)$

Cuando la media es un múltiplo de $\theta=(0, ..., 0, 1)$. Disminuímos la dimensión para que se aprecie más.

Vayamos ahora al caso donde la matriz de covarianzas es un múltiplo desconocido de la identidad. Tendremos que aproximarla mediante una distribució ji cuadrado.

Comparando con otras distribuciones

t-Student

Comparamos ahora el estimador de James-Stein entre el caso normal i el t-Student, suponiendo el vector de medias $\theta=0$.

Ahora hacemos lo mismo pero en el caso en que la norma del vector de medias aumenta. Establecemos 6 grados de libertad, por ejemplo.

Veamos si también tenemos una contración por ejemplo con distribuciones de Poisson. Lo comparamos al caso Gaussiano

Poisson

Ji cuadrado

Vemos que el estimador de James-Stein también produce cierta mejora. Fijamos la dimensión $p$ y establecemos por ejemplo $6$ grados de libertad:

Y en función de la dimensión $p$:

Vemos también que fijando la dimensión $p$, la mejora se mantiene para distintos grados de libertad: