¿Cuál sería un modelo bayesiano robusto para estimar la escala de una distribución más o menos normal?

32

Existe una serie de estimadores robustos de escala . Un ejemplo notable es la mediana de la desviación absoluta que se relaciona con la desviación estándar como σ=METROUNAre1.4826 . En un marco bayesiano, existen varias formas de estimar de manera sólida la ubicación de una distribución más o menos normal (por ejemplo, una Normal contaminada por valores atípicos), por ejemplo, se podría suponer que los datos se distribuyen como en distribución o distribución de Laplace. Ahora mi pregunta:

¿Cuál sería un modelo bayesiano para medir la escala de una distribución más o menos normal de una manera robusta, robusta en el mismo sentido que el MAD o estimadores robustos similares?

Como es el caso con MAD, sería genial si el modelo bayesiano pudiera acercarse al SD de una distribución normal en el caso en que la distribución de los datos en realidad se distribuye normalmente.

editar 1:

Un ejemplo típico de un modelo que es robusto contra la contaminación / valores atípicos al asumir los datos yyo es más o menos normal es usar en la distribución como:

yit(m,s,ν)

Donde m es la media, s es la escala y ν es el grado de libertad. Con priors adecuados en m,s y ν , m será una estimación de la media de yi que va a ser robusto frente a valores atípicos. Sin embargo, s no será una estimación consistente de la SD de yi ya que s depende de ν . Por ejemplo, si ν se fijaría en 4.0 y el modelo anterior se ajustaría a una gran cantidad de muestras de un distribución entonces s estaría alrededor de 0.82. Lo que estoy buscando es un modelo que sea robusto, como el modelo t, pero para el SD en lugar de (o además de) la media.Norm(μ=0,σ=1)s

editar 2:

Aquí sigue un ejemplo codificado en R y JAGS de cómo el modelo t 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 
Rasmus Bååth
fuente
Tal vez no sea lo suficientemente robusto, pero la distribución de chi-cuadrado es el conjugado elegido generalmente antes de la inversa de la varianza.
Mike Dunlavey
Es posible que desee ver si la primera respuesta a esta pregunta stats.stackexchange.com/questions/6493/… es suficiente para usted; puede que no sea así, pero tal vez lo sea.
jbowman
¿Cuál es su nivel de contaminación previo? ¿La contaminación será sistemática? ¿Aleatorio? ¿Será generado por una sola distribución o múltiples distribuciones? ¿Sabemos algo sobre la distribución del ruido? Si se conocen al menos algunas de las cosas anteriores, entonces podríamos adaptarnos a algún tipo de modelo de mezcla. De lo contrario, no estoy seguro de cuáles son realmente sus creencias sobre este problema, y ​​si no tiene ninguna, esto parece un entorno muy vago. Debe arreglar algo, de lo contrario puede elegir un punto al azar y declararlo como el único punto generado con Gauss.
significado para el
Pero, en general, podría ajustarse a una distribución t que sea más resistente contra los valores atípicos, o una mezcla de distribuciones t. Estoy seguro de que hay muchos documentos, aquí hay uno de Bishop research.microsoft.com/en-us/um/people/cmbishop/downloads/… y aquí hay un paquete R para las mezclas: maths.uq.edu. au / ~ gjm / mix_soft / EMMIX_R / EMMIX-manual.pdf
significa significado el
1
Su es cierto para una población normalmente distribuida, pero no para la mayoría de las otras distribucionesσ=MAD1.4826
Henry

Respuestas:

10

La inferencia bayesiana en un modelo de ruido T con un previo apropiado dará una estimación sólida de ubicación y escala. Las condiciones precisas que la probabilidad y la necesidad previa de satisfacer se dan en el documento modelo de robustez bayesiano de los parámetros de ubicación y escala de Andrade y O'Hagan (2011). Las estimaciones son sólidas en el sentido de que una sola observación no puede hacer que las estimaciones sean arbitrariamente grandes, como se demuestra en la figura 2 del documento.

Cuando los datos se distribuyen normalmente, la SD de la distribución T ajustada (para fijo ) no coincide con la SD de la distribución generadora. Pero esto es fácil de arreglar. Sea σ la desviación estándar de la distribución generadora y s sea ​​la desviación estándar de la distribución T ajustada. Si los datos se escalan en 2, entonces, por la forma de la probabilidad, sabemos que s debe escalar en 2. Esto implica que s = σ f ( ν ) para alguna función fija f . Esta función se puede calcular numéricamente mediante simulación a partir de una normal estándar. Aquí está el código para hacer esto:νσsss=σf(ν)f

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 obtengo f ( ν ) = 1.18 . El estimador deseado es entonces σ = s / f ( ν ) .ν=4f(ν)=1.18σ^=s/f(ν)

Tom Minka
fuente
1
Buena respuesta (+1). "en el sentido de que una sola observación no puede hacer que las estimaciones sean arbitrariamente grandes", por lo que el punto de ruptura es 2 / n (me preguntaba sobre esto) ... Como punto de comparación, para el procedimiento ilustrado en mi respuesta es n / 2.
user603
¡Wow gracias! Pregunta de seguimiento difusa. ¿Tendría sentido "corregir" la escala para que sea coherente con la SD en el caso Normal? El caso de uso en el que estoy pensando es al informar una medida de propagación. No tendría ningún problema con la escala de informes, pero sería bueno informar algo que sería coherente con SD ya que es la medida más común de propagación (al menos en psicología). ¿Ves una situación en la que esta corrección conduciría a estimaciones extrañas e inconsistentes?
Rasmus Bååth
6

Cuando haga una pregunta sobre un problema muy preciso (estimación sólida), le ofreceré una respuesta igualmente precisa. Primero, sin embargo, comenzaré a tratar de disipar una suposición injustificada. No es cierto que haya una estimación bayesiana robusta de la ubicación (hay estimadores bayesianos de ubicaciones, pero como ilustraré a continuación no son robustos y, aparentemente , incluso el estimador robusto de ubicación más simple no es bayesiano). En mi opinión, las razones de la ausencia de superposición entre el paradigma 'bayesiano' y el 'robusto' en el caso de la ubicación explica en gran medida por qué tampoco hay estimadores de dispersión que sean robustos y bayesianos.

Con priors adecuados en y ν , m será una estimación de la media de y im,sνmyi que va a ser robusto frente a valores atípicos.

En realidad no. Las estimaciones resultantes solo serán sólidas en un sentido muy débil de la palabra robusto. Sin embargo, cuando decimos que la mediana es robusta para los valores atípicos, nos referimos a la palabra robusta en un sentido mucho más fuerte. Es decir, en estadísticas sólidas, la solidez de la mediana se refiere a la propiedad de que si calcula la mediana en un conjunto de datos de observaciones extraídas de un modelo continuo unimodal y luego reemplaza menos de la mitad de estas observaciones por valores arbitrarios , el valor de la mediana calculada en los datos contaminados está cerca del valor que hubiera tenido si lo hubiera calculado en el conjunto de datos original (no contaminado). Entonces, es fácil demostrar que la estrategia de estimación que propone en el párrafo que cité anteriormente definitivamente no es robusto en el sentido de cómo la palabra se entiende típicamente para la mediana.

No estoy completamente familiarizado con el análisis bayesiano. Sin embargo, me preguntaba qué hay de malo en la siguiente estrategia, ya que parece simple, eficaz y, sin embargo, no se ha considerado en las otras respuestas. Lo anterior es que la buena parte de los datos se extrae de una distribución simétrica y que la tasa de contaminación es inferior a la mitad. Entonces, una estrategia simple sería:F

  1. Calcule la mediana / locura de su conjunto de datos. Luego calcule:
    zi=|ximed(x)|mad(x)
  2. excluya las observaciones para las cuales (este es el cuantil α de la distribución de z cuando x F ). Esta cantidad está disponible para muchas opciones de F y puede arrancarse para los demás.zi>qα(z|xF)αzxFF
  3. Ejecute un análisis bayesiano (habitual, no robusto) sobre las observaciones no rechazadas.

EDITAR:

Gracias al OP por proporcionar un código R autónomo para realizar un análisis bayesiano bonna fide del problema.

El siguiente código compara el enfoque bayesiano sugerido por el OP con su alternativa de la literatura estadística robusta (por ejemplo, el método de ajuste propuesto por Gauss para el caso donde los datos pueden contener tanto como n/22 valores atípicos y la distribución de buena parte de los datos es gaussiana).

La parte central de los datos es :N(1000,1)

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

Agregue 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 el valor 1 para los valores atípicos. Comienzo con el enfoque sugerido por el OP:

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)

Yo 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 

(silencioso lejos de los valores objetivo)

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 objetivo)

zthF
t

  • [1] Ricardo A. Maronna, Douglas R. Martin, Victor J. Yohai (2006). Estadísticas robustas: teoría y métodos (Serie Wiley en probabilidad y estadística).
  • Huber, PJ (1981). Estadísticas robustas. Nueva York: John Wiley and Sons.
usuario603
fuente
1
Bueno, la t a menudo se propone como una alternativa robusta a la distribución normal. No sé si esto es en sentido débil o no. Ver por ejemplo: Lange, KL, Little, RJ y Taylor, JM (1989). Modelado estadístico robusto utilizando la distribución t. Revista de la Asociación Americana de Estadística , 84 (408), 881-896. pdf
Rasmus Bååth
1
Este es el sentido débil. Si tiene un código R que implementa el procedimiento que sugiere, me complacerá ilustrar mi respuesta con un ejemplo. de lo contrario, puede obtener más explicaciones en el capítulo 2 de este libro de texto.
user603
El procedimiento que sugiero se describe básicamente aquí: indiana.edu/~kruschke/BEST, incluido el código R. ¡Tendré que pensar en su solución! Sin embargo, no parece Bayesiano en el sentido de que no modela todos los datos, solo el subconjunto que "sobrevive" al paso 2.
Rasmus Bååth
¡Gracias por tu interesante discusión! Sin embargo, su respuesta no es lo que busco, porque (1) no describe un procedimiento bayesiano, describe más un paso de preparación de datos sobre cómo eliminar los valores atípicos (2) su procedimiento no resulta en un estimador consistente de DE , es decir, si toma muestras de una distribución normal y el número de puntos de datosno se acercará al SD "verdadero", sino que su estimación será un poco baja. Tampoco compro completamente su definición de robusto (su definición no es como la he visto en la mayoría de la literatura bayesiana que he encontrado)
Rasmus Bååth
1
¡Ya lo hice!
Rasmus Bååth
1

En el análisis bayesiano, el uso de la distribución inversa de Gamma como prioridad para la precisión (el inverso de la varianza) es una opción común. O la distribución inversa de Wishart para modelos multivariados. Agregar un previo a la variación mejora la robustez frente a los valores atípicos.

Hay un buen artículo de Andrew Gelman: "Distribuciones previas para parámetros de varianza en modelos jerárquicos" en el que discute qué buenas elecciones pueden ser para los previos en las varianzas.

jpmuc
fuente
44
Lo siento, pero no veo cómo esto responde la pregunta. No pedí un previo robusto, sino un modelo robusto .
Rasmus Bååth
0

Un estimador robusto para el parámetro de ubicación μ de algún conjunto de datos de tamaño norte se obtiene cuando se asigna un Jeffreys antes de la varianza σ2 de la distribución normal, y calcula el marginal para μ, produciendo un t distribución con norte grados de libertad.

Del mismo modo, si desea un estimador robusto para la desviación estándar σ de algunos datos re, podemos hacer lo siguiente:

Primero, suponemos que los datos se distribuyen normalmente cuando se conocen su media y la desviación estándar. Por lo tanto,

reEl |μ,σN(μ,σ2)
and if D(d1,,dN) then
p(D|μ,σ2)=1(2πσ)Nexp(N2σ2((mμ2)+s2))
where the sufficient statistics m and s2 are
m=1Ni=1Ndis2=1Ni=1Ndi2m2
In addition, using Bayes' theorem, we have
p(μ,σ2|D)p(D|μ,σ2)p(μ,σ2)
A convenient prior for (μ,σ2) is the Normal-invese-gamma family, which covers a wide range of shapes and is conjugate to this likelihood. This means that the posterior distribution p(μ,σ2|D) still belongs to the normal-inverse-gamma family, and its marginal p(σ2|D) is an inverse gamma distribution parameterized as
σ2|DIG(α+N/2,2β+Ns2)α,β>0
From this distribution, we can take the mode, which will give us an estimator for σ2. This estimator will be more or less tolerant to small excursions from misspecifications on the model by varying α and/or β. The variance of this distribution will then provide some indication on the fault-tolerance of the estimate. Since the tails of the inverse gamma are semi-heavy, you get the kind of behaviour you would expect from the t distribution estimate for μ that you mention.
yannick
fuente
1
"A robust estimator for the location parameter μ of some dataset of size N is obtained when one assigns a Jeffreys prior to the variance σ2 of the normal distribution." Isn't this Normal model you describe a typical example of a non-robust model? That is, a single value that is off can have great influence on the parameters of the model. There is a big difference between the posterior over the mean being a t-distribution (as in your case) and the distribution for the data being a t-distribution (as is a common example of a robust Bayesian model for estimating the mean).
Rasmus Bååth
1
It all depends on what you mean by robust. What you are saying right now is that you would like robustness wrt data. What I was proposing was robustness wrt model mis-specification. They are both different types of robustness.
yannick
2
I would say that the examples I gave, MAD and using a t distribution as the distribution for the data are examples of robustness with respect to data.
Rasmus Bååth
I would say Rasmus is right and so would Gelman er al in BDA3, as would a basic understanding that th t distribution has fatter tails than the normal for the same location parameter
Brash Equilibrium
0

I have followed the discussion from the original question. Rasmus when you say robustness I am sure you mean in the data (outliers, not miss-specification of distributions). I will take the distribution of the data to be Laplace distribution instead of a t-distribution, then as in normal regression where we model the mean, here we will model the median (very robust) aka median regression (we all know). Let the model be:

Y=βX+ϵ, ϵ has laplace(0,σ2).

Of course our goal is to estimate model parameters. We expect our priors to be vague to have an objective model. The model at hand has a posterior of the form f(β,σ,Y,X). Giving β a normal prior with large variance makes such a prior vague and a chis-squared prior with small degrees of freedom to mimic a jeffrey's prior(vague prior) is given to to σ2. With a Gibbs sampler what happens? normal prior+laplace likehood=???? we do know. Also chi-square prior +laplace likelihood=??? we do not know the distribution. Fortunately for us there is a theorem in (Aslan,2010) that transforms a laplace likelihood to a scale mixture of normal distributions which then enable us to enjoy the conjugate properties of our priors. I think the whole process described is fully robust in terms of outliers. In a multivariate setting chi-square becomes a a wishart distribution, and we use multivariate laplace and normal distributions.

Chamberlain Foncha
fuente
2
Your solution seems to be focused on robust estimation of the location(mean/median). My question was rather about estimation of scale with the property of consistency with respect to retrieving the SD when the data generating distribution actually is normal.
Rasmus Bååth
With a robust estimate of the location, the scale as function of the location immediately benefits from the robustness of the location. There is no other way of making the scale robust.
Chamberlain Foncha
Anyway I must say I am eagerly waiting to see how this problem will be tackled most especially with a normal distribution as you emphasized.
Chamberlain Foncha
0

Suppose that you have K groups and you want to model the distribution of their sample variances, perhaps in relation to some covariates x. That is, suppose that your data point for group k1K is Var(yk)[0,). The question here is, "What is a robust model for the likelihood of the sample variance?" One way to approach this is to model the transformed data ln[Var(yk)] as coming from a t distribution, which as you have already mentioned is a robust version of the normal distribution. If you don't feel like assuming that the transformed variance is approximately normal as n, then you could choose a probability distribution with positive real support that is known to have heavy tails compared to another distribution with the same location. For example, there is a recent answer to a question on Cross Validated about whether the lognormal or gamma distribution has heavier tails, and it turns out that the lognormal distribution does (thanks to @Glen_b for that contribution). In addition, you could explore the half-Cauchy family.

Similar reasoning applies if instead you are assigning a prior distribution over a scale parameter for a normal distribution. Tangentially, the lognormal and inverse-gamma distributions are not advisable if you want to form a boundary avoiding prior for the purposes of posterior mode approximation because they peak sharply if you parameterize them so that the mode is near zero. See BDA3 chapter 13 for discussion. So in addition to identifying a robust model in terms of tail thickness, keep in mind that kurtosis may matter to your inference, too.

I hope this helps you as much as your answer to one of my recent questions helped me.

Brash Equilibrium
fuente
1
My question was about the situation when you have one group and how to robustly estimate the scale of that group. In the case of outliers I don't believe the sample variance is considered robust.
Rasmus Bååth
If you have one group, and you are estimating its normal distribution, then your question applies to the form of the prior over its scale parameter. As my answer implies, you can use a t distribution over its log transformation or choose a fat tailed distribution with positive real support, being careful about other aspects of that distribution such as its kurtosis. Bottom line, if you wan a robust model for a scale parameter, use a t distribution over its log transform or some other fat tailed distribution.
Brash Equilibrium