¿Qué sería de un robusto modelo Bayesiano para estimar la escala de aproximadamente una distribución normal?

Question

Existe una serie de estimadores robustos de escala. Un ejemplo notable es la mediana de la desviación absoluta que se refiere a la desviación estándar como $\sigma = \mathrm{MAD}\cdot1.4826$. En un marco Bayesiano existen un número de maneras de firmeza a la estimación de la ubicación de aproximadamente una distribución normal (digamos Normal contaminados por los valores atípicos), por ejemplo, se podría suponer que los datos se distribuye como una t de la distribución o distribución de Laplace. Ahora mi pregunta:

¿Qué sería de un modelo Bayesiano para la medición de la escala de aproximadamente una distribución normal en una forma robusta de ser, robusto en el mismo sentido como el LOCO o similar estimadores robustos?

Como es el caso de la LOCA, que sería genial si el modelo Bayesiano podría acercarse a la SD de una distribución normal, en el caso de la distribución de los datos está normalmente distribuida.

edit 1:

Un ejemplo típico de un modelo que es robusto contra la contaminación/valores atípicos cuando asumiendo que los datos $y_i$ es aproximadamente normal es el uso de una distribución t como:

$$y_i \sim \mathrm{t}(m, s,\nu)$$

Donde $m$ es decir, $s$ es la escala, y $\nu$ es el grado de libertad. Con adecuado de los priores en $m, s$ y $\nu$, $m$ será una estimación de la media de $y_i$ que va a ser robusto frente a los valores atípicos. Sin embargo, $s$ no será una estimación consistente de la SD de $y_i$ $s$ depende de $\nu$. Por ejemplo, si $\nu$ sería fijo a 4.0 y el modelo de arriba podría ser instalado en un gran número de muestras de un $\mathrm{Norm}(\mu=0,\sigma=1)$ de distribución, a continuación, $s$ sería de alrededor de 0.82. Lo que estoy buscando es un modelo robusto, como el modelo t, pero para la SD en lugar de (o además de) de la media.

edit 2:

Aquí sigue un código de ejemplo en R y PUNTAS de cómo el t-modelo mencionado anteriormente es más robusto con respecto a la media.

# generating some contaminated data
y <- c( rnorm(100, mean=10, sd=10), 
        rnorm(10, mean=100, sd= 100))

#### A "standard" normal model ####
model_string <- "model{
  for(i in 1:length(y)) {
    y[i] ~ dnorm(mu, inv_sigma2)
  }

  mu ~ dnorm(0, 0.00001)
  inv_sigma2 ~ dgamma(0.0001, 0.0001)
  sigma <- 1 / sqrt(inv_sigma2)
}"

model <- jags.model(textConnection(model_string), list(y = y))
mcmc_samples <- coda.samples(model, "mu", n.iter=10000)
summary(mcmc_samples)

### The quantiles of the posterior of mu
##  2.5%   25%   50%   75% 97.5% 
##   9.8  14.3  16.8  19.2  24.1 

#### A (more) robust t-model ####
library(rjags)
model_string <- "model{
  for(i in 1:length(y)) {
    y[i] ~ dt(mu, inv_s2, nu)
  }

  mu ~ dnorm(0, 0.00001)
  inv_s2 ~ dgamma(0.0001,0.0001)
  s <- 1 / sqrt(inv_s2)
  nu ~ dexp(1/30) 
}"

model <- jags.model(textConnection(model_string), list(y = y))
mcmc_samples <- coda.samples(model, "mu", n.iter=1000)
summary(mcmc_samples)

### The quantiles of the posterior of mu
## 2.5%   25%   50%   75% 97.5% 
##8.03  9.35  9.99 10.71 12.14 

Answer 1

Inferencia bayesiana en un T de ruido, con un modelo apropiado antes, se le dará una estimación robusta de localización y escala. Las condiciones precisas que la probabilidad y antes de que se deben satisfacer son dadas en el papel robustez Bayesiana de la modelización de la ubicación y la escala de los parámetros por Andrade y O'Hagan (2011). Las estimaciones son robustas en el sentido de que una sola observación no se pueden hacer las estimaciones arbitrariamente grandes, como se muestra en la figura 2 de la hoja de papel.

Cuando los datos están distribuidos normalmente, la tarjeta SD de la equipada con T de distribución (fijo $\nu$) no coincide con la SD de la generación de la distribución. Pero esto es fácil de arreglar. Deje que $\sigma$ ser la desviación estándar de la generación de la distribución y deje de $s$ ser la desviación estándar de los armarios T de distribución. Si los datos se ajustan a 2, luego de la forma de la probabilidad sabemos que $$ s debe escala 2. Esto implica que $s = \sigma f(\nu)$ para algunos de la función fija de $f$. Esta función puede ser calculada numéricamente mediante la simulación de una normal estándar. Aquí está el código para hacer esto:

library(stats)
library(stats4)
y = rnorm(100000, mean=0,sd=1)
nu = 4
nLL = function(s) -sum(stats::dt(y/s,nu,log=TRUE)-log(s))
fit = mle(nLL, start=list(s=1), method="Brent", lower=0.5, upper=2)
# the variance of a standard T is nu/(nu-2)
print(coef(fit)*sqrt(nu/(nu-2)))

Por ejemplo, en $\nu=4$ I get $f(\nu)=1.18$. El deseado estimador es entonces $\hat{\sigma} = s/f(\nu)$.

Answer 2

Como usted está haciendo una pregunta sobre un problema muy concreto (estimación robusta), te ofrecen igualmente una respuesta precisa. En primer lugar, sin embargo, voy a empezar a estar tratando de disipar un injustificado de la asunción. No es cierto que existe una fuerte bayesiano estimación de la ubicación (hay estimadores bayesianos de lugares, pero como se ilustra más abajo no son robustos y, al parecer, incluso el más simple estimador robusto de la ubicación no está bayesiano) . En mi opinión, las razones de la ausencia de solapamiento entre la 'bayesiano' y 'robusto' paradigma en el caso de la posición va un largo camino en la explicación de por qué no también hay estimadores de dispersión que son a la vez robusto y bayesiano .

Con adecuado de los priores en $m, s$ y $\nu$, $m$ será una estimación de la media de $y_i$ que va a ser robusto frente a los valores atípicos.

En realidad, no. Las estimaciones resultantes sólo será robusta en un muy débil sentido de la palabra robusto. Sin embargo, cuando decimos que la mediana es robusto a los valores atípicos, nos referimos a la palabra robusto en mucho más sentido. Que es, en estadísticas robustas, la solidez de la mediana se refiere a la propiedad de que si usted calcular la mediana en un conjunto de datos de observaciones extraídas de una uni-modal, modelo continuo y, a continuación, reemplace menos de la mitad de estas observaciones arbitraria de valores, el valor de la mediana calculada sobre la contaminación de los datos es cercano al valor que hubiera tenido había calculado que en el original (no contaminado) conjunto de datos. Entonces, es fácil demostrar que la estrategia de estimación que usted propone en el párrafo que he citado anteriormente es definitivamente no es robusto en el sentido de cómo la palabra generalmente se entiende por la mediana.

Estoy totalmente familiarizado con el análisis Bayesiano. Sin embargo, me preguntaba ¿qué está mal con la siguiente estrategia como parece simple, eficaz y, sin embargo, no se ha considerado en las otras respuestas. La anterior es que la buena parte de los datos provienen de una distribución simétrica de $F$ y que la tasa de contaminación es menos de la mitad. Entonces, una estrategia simple sería:

1. calcular la mediana/mad de su conjunto de datos. A continuación, calcular: $$z_i=\frac{|x_i-\mbox{med}(x)|}{\mbox{mad}(x)}$$
2. excluir las observaciones de los cuales $z_i>q_{\alpha}(z|x\sim F)$ (este es el $\alpha$ cuantil de la distribución de los $z$ cuando $x\sim F$). Esta cantidad está disponible para muchos la elección de $F$ y puede ser bootstrap para los demás.
3. Ejecutar una (como es habitual, no sólido) análisis Bayesiano sobre la no-rechazó las observaciones.

EDITAR:

Gracias a la OP para la prestación de un auto contenida R código de conducta de un bonna fide análisis bayesiano de la problema.

el código siguiente se compara el enfoque bayesiano sugerido por la O. P. a es la alternativa de las estadísticas robustas de la literatura (por ejemplo, el método de ajuste propuesto por Gauss para el caso en que los datos pueden contener tanto como $n/2-2$ los valores atípicos y la distribución de la parte buena de los datos es de Gauss).

parte central de los datos es de $\mathcal{N}(1000,1)$:

n<-100
set.seed(123)
y<-rnorm(n,1000,1)

Añadir una cierta cantidad de contaminantes:

y[1:30]<-y[1:30]/100-1000 
w<-rep(0,n)
w[1:30]<-1

el índice w toma valor 1 para los valores atípicos. Empiezo con el enfoque sugerido por la O. P.:

library("rjags")
model_string<-"model{
  for(i in 1:length(y)){
    y[i]~dt(mu,inv_s2,nu)
  }
  mu~dnorm(0,0.00001)
  inv_s2~dgamma(0.0001,0.0001)
  s<-1/sqrt(inv_s2)
  nu~dexp(1/30) 
}"

model<-jags.model(textConnection(model_string),list(y=y))
mcmc_samples<-coda.samples(model,"mu",n.iter=1000)
print(summary(mcmc_samples)$statistics[1:2])
summary(mcmc_samples)

Obtengo:

     Mean        SD 
384.2283  97.0445 

y:

2. Quantiles for each variable:

 2.5%   25%   50%   75% 97.5% 
184.6 324.3 384.7 448.4 577.7 

(tranquilo, lejos por tanto de los valores de destino)

Para el método robusto,

z<-abs(y-median(y))/mad(y)
th<-max(abs(rnorm(length(y))))
print(c(mean(y[which(z<=th)]),sd(y[which(z<=th)])))

uno obtiene:

 1000.149 0.8827613

(muy cerca de los valores de destino)

El segundo resultado es mucho más cerca de los valores reales. Pero se pone peor. Si se clasifican como valores atípicos esas observaciones para que el estimado de $z$-score es mayor que th (recordemos que el antes es que $F$ es Gaussiano), a continuación, el enfoque bayesiano se encuentra que todas las observaciones son valores atípicos (la robusta procedimiento, en contraste, las banderas de todos y sólo los valores extremos como tal). Esto también implica que si se ejecuta una costumbre (no-robusto) análisis bayesiano en los datos no se clasifican como valores atípicos por la robusta procedimiento, se debe hacer bien (por ejemplo, cumplir con los objetivos establecidos en su pregunta).
Este es sólo un ejemplo, pero en realidad es bastante sencillo demostrar que (y no se puede hacer formalmente, véase, por ejemplo, en el capítulo 2 de [1]) los parámetros de un estudiante $t$ distribución de módulos a los datos que no se puede depender para revelar los valores atípicos.

• [1]Ricardo A. Maronna, Douglas R. Martin, Victor J. Yohai (2006). Estadísticas robustas: Teoría y Métodos (Wiley Serie en Probabilidad y Estadística).
• Huber, P. J. (1981). Estadísticas Robustas. Nueva York: John Wiley and Sons.

Answer 3

En el análisis bayesiano utilizando la inversa de la distribución Gamma como antes de la precisión (la inversa de la varianza) es una opción común. O la inversa de la distribución de Wishart para los modelos multivariables. La adición de un previo en la variación mejora la robustez en contra de los valores atípicos.

Hay un bonito papel de Andrew Gelman: "Antes de las distribuciones de variación de parámetros en los modelos jerárquicos" donde se discute qué buenas opciones para los priores de las desviaciones pueden ser.

Answer 4

Un robusto estimador del parámetro de localización de $\mu$ de algún conjunto de datos de tamaño $N$ se obtiene cuando se asigna un Jeffreys antes de la varianza de $\sigma^2$ de la distribución normal, y se calcula la marginal de $\mu$, produciendo un $t$ de distribución con $$ N grados de libertad.

Del mismo modo, si usted quiere un robusto estimador de la desviación estándar de $\sigma$ de algunos datos $D$, podemos hacer lo siguiente:

En primer lugar, suponemos que los datos están distribuidos normalmente, cuando su media y desviación estándar son conocidos. Por lo tanto, $$\left.D\right|_{\mu\sigma} \sim \mathcal{N}(\mu\sigma^2)$$ y si $D \equiv (d_1,\ldots,d_N)$ $$p(D|\mu\sigma^2) = \frac{1}{(\sqrt{2\pi}\sigma)^N} \exp\left(-\frac{N}{2\sigma^2}\left((m-\mu^2)+s^2\right)\right)$$ donde la suficiente estadísticas $m$ y $s^2$ son $$m=\frac{1}{N}\sum_{i=1}^N d_i \quad s^2 = \frac{1}{N}\sum_{i=1}^N d_i^2 m^2$$ Además, usando el teorema de Bayes, tenemos $$p(\mu\sigma^2|D) \propto p(D|\mu\sigma^2) p(\mu\sigma^2)$$ Un conveniente antes de $(\mu\sigma^2)$ es la Normal-invese-gamma de la familia, que cubre una amplia gama de formas y es conjugado a esta probabilidad. Esto significa que la distribución posterior $p(\mu\sigma^2|D)$ aún pertenece a la normal inverso-gamma de la familia, y su marginales $p(\sigma^2|D)$ es una inversa de la distribución gamma con parámetros como $$\left.\sigma^2\right|_{D} \sim \mathcal{IG}\left(\alpha+N/2,2\beta+Ns^2\right) \qquad \alpha,\beta>0$$ A partir de esta distribución, se puede tomar el modo, que nos dará un estimador para $\sigma^2$. Este estimador será más o menos tolerantes a las pequeñas excursiones de misspecifications en el modelo mediante la variación de $\alpha$ y/o $\beta$. La varianza de esta distribución, a continuación, proporcionar alguna indicación de la tolerancia a fallos de la estimación. Desde las colas de la inversa de la gamma son semi-pesados, se obtiene el tipo de comportamiento que se espera de el $t$ distribución de la estimación para $\mu$ que usted menciona.

Answer 5

He seguido la discusión de la pregunta original. Rasmus cuando dices robustez estoy seguro de que la media de los datos (valores atípicos, no se pierda la especificación de las distribuciones). Voy a tomar la distribución de los datos a ser de Laplace de distribución en lugar de una distribución t, entonces como en condiciones normales de regresión donde el modelo de la media, aquí vamos a modelar la mediana (muy sólida) también conocido como la mediana de la regresión (todos sabemos). Deje que el modelo de ser:

Y=\beta X+\epsilon$,$\epsilon$ha laplace$(0,$$\sigma^2)$.

Por supuesto, nuestro objetivo es estimar los parámetros del modelo. Esperamos que nuestros priores que ser vago para tener un modelo objetivo. El modelo en la mano tiene una parte posterior de la forma $f(\beta\sigma,Y,X)$. Dando $\beta$ previa normal con varianza grande hace que sea un antes vago y un chis-squared antes con pequeños grados de libertad para imitar una jeffrey previo(vago antes de) a $\sigma^2$. Con un muestreador de Gibbs, ¿qué sucede? normal antes+de laplace likehood=???? hacemos saber. También de chi-cuadrado antes +de laplace, la probabilidad=??? no sabemos la distribución. Afortunadamente para nosotros no es un teorema (Aslan,2010) que transforma una de laplace, la probabilidad de una escala de mezcla de distribuciones normales que luego nos permiten disfrutar de las propiedades del conjugado de nuestras prioridades. Creo que todo el proceso descrito es totalmente robusta en términos de valores atípicos. En un multivariante de ajuste de chi-cuadrado se convierte en un una distribución de wishart, y utilizamos multivariante de laplace y distribuciones normales.