Existe um número de estimadores robustos de escala . Um exemplo notável é o desvio absoluto mediana que se relaciona com o desvio padrão como $\sigma = \mathrm{MAD}\cdot1.4826$ . Em uma estrutura bayesiana, existem várias maneiras de estimar com robustez a localização de uma distribuição aproximadamente normal (por exemplo, um Normal contaminado por outliers); por exemplo, pode-se supor que os dados sejam distribuídos na distribuição ou na distribuição de Laplace. Agora minha pergunta:

O que seria um modelo bayesiano para medir a escala de uma distribuição aproximadamente normal de uma maneira robusta, robusta no mesmo sentido que o MAD ou estimadores robustos similares?

Como é o caso do MAD, seria interessante se o modelo bayesiano pudesse abordar o DP de uma distribuição normal no caso em que a distribuição dos dados realmente é normalmente distribuída.

editar 1:

Um exemplo típico de um modelo que é robusto contra contaminação / outliers ao assumir que os dados são aproximadamente normais está usando na distribuição como:$y_i$

Onde é a média, é a escala e é o grau de liberdade. Com antecedentes adequados em e , será uma estimativa da média de que vai ser robusta contra os outliers. No entanto, não será uma estimativa consistente do DP de pois depende de . Por exemplo, se fosse fixado em 4,0 e o modelo acima fosse ajustado a um grande número de amostras de uma distribuição então$m$$s$$\nu$$m, s$$\nu$$m$$y_i$$s$$y_i$$s$$\nu$$\nu$$\mathrm{Norm}(\mu=0,\sigma=1)$$s$seria em torno de 0,82. O que estou procurando é um modelo robusto, como o modelo t, mas para o SD em vez de (ou além) da média.

editar 2:

Segue um exemplo codificado em R e JAGS de como o modelo t mencionado acima é mais robusto em relação à média.

# 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 

A inferência bayesiana em um modelo de ruído T com um prévio apropriado fornecerá uma estimativa robusta de localização e escala. As condições precisas que a probabilidade e a necessidade prévia precisam satisfazer são fornecidas no artigo Modelagem da robustez bayesiana dos parâmetros de localização e escala por Andrade e O'Hagan (2011). As estimativas são robustas no sentido de que uma única observação não pode tornar as estimativas arbitrariamente grandes, como demonstrado na figura 2 do artigo.

Quando os dados são normalmente distribuídos, o SD da distribuição T ajustada (para fixo ) não corresponde ao SD da distribuição geradora. Mas isso é fácil de consertar. Deixe σ ser o desvio-padrão da distribuição de geração e deixar de ser o desvio-padrão da distribuição t equipada. Se os dados forem redimensionados por 2, então, a partir da forma da probabilidade, sabemos que s deve ser redimensionado por 2. Isso implica que s = σ f ( v ) para alguma função fixa f . Esta função pode ser calculada numericamente por simulação a partir de um normal padrão. Aqui está o código para fazer isso:$\nu$$\sigma$$s$$s$$s = \sigma f(\nu)$$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 exemplo, em , obtenho f ( ν ) = 1,18 . O estimador desejado é, em seguida, σ = s / f ( ν ) .$\nu=4$$f(\nu)=1.18$$\hat{\sigma} = s/f(\nu)$

Como você está fazendo uma pergunta sobre um problema muito preciso (estimativa robusta), vou oferecer uma resposta igualmente precisa. Primeiro, no entanto, começarei a tentar dissipar uma suposição injustificada. Não é verdade que exista uma estimativa bayesiana robusta de localização (existem estimadores bayesianos de localizações, mas, como ilustro abaixo, eles não são robustos e, aparentemente , mesmo o estimador robusto mais simples de localização não é bayesiano). Na minha opinião, as razões para a ausência de sobreposição entre o paradigma 'bayesiano' e 'robusto' no caso de localização explicam por que não existem estimadores de dispersão robustos e bayesianos.

Com antecedentes adequados em e ν , m será uma estimativa da média de y i que vai ser robusta contra os outliers.$m, s$$\nu$$m$$y_i$

Na verdade não. As estimativas resultantes serão robustas apenas em um sentido muito fraco da palavra robusto. No entanto, quando dizemos que a mediana é robusta para valores extremos, queremos dizer a palavra robusto em um sentido muito mais forte. Ou seja, em estatísticas robustas, a robustez da mediana refere-se à propriedade que, se você calcular a mediana em um conjunto de dados de observações extraídas de um modelo contínuo uni-modal e depois substituir menos da metade dessas observações por valores arbitrários , o valor da mediana calculada nos dados contaminados é próximo ao valor que você teria se o tivesse calculado no conjunto de dados original (não contaminado). Então, é fácil mostrar que a estratégia de estimativa que você propõe no parágrafo citado acima definitivamente não é robusto no sentido de como a palavra é tipicamente entendida para a mediana.

Não estou familiarizado com a análise bayesiana. No entanto, eu queria saber o que há de errado com a estratégia a seguir, pois parece simples, eficaz e ainda não foi considerada nas outras respostas. O anterior é que a boa parte dos dados é extraída de uma distribuição simétrica e que a taxa de contaminação é menor que a metade. Então, uma estratégia simples seria:$F$

1. calcule a mediana / louco do seu conjunto de dados. Depois calcule: $$z_i=\frac{|x_i-\mbox{med}(x)|}{\mbox{mad}(x)}$$
2. excluir as observações de que (esta é a α quantil da distribuição de Z quando X ~ F ). Essa quantidade está disponível para muitas opções de F e pode ser inicializada para as outras.$z_i>q_{\alpha}(z|x\sim F)$$\alpha$$z$$x\sim F$$F$
3. Execute uma análise bayesiana (usual, não robusta) nas observações não rejeitadas.

EDITAR:

Agradecimentos ao OP por fornecer um código R independente para conduzir uma análise bayesiana de boa-fé do problema.

o código abaixo compara a abordagem bayesiana sugerida pelo OP à sua alternativa da literatura estatística robusta (por exemplo, o método de ajuste proposto por Gauss para o caso em que os dados podem conter até outliers e a distribuição do boa parte dos dados é gaussiana).$n/2-2$

a parte central dos dados é :$\mathcal{N}(1000,1)$

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

Adicione uma certa quantidade de contaminantes:

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

o índice w assume o valor 1 para os valores discrepantes. Começo com a abordagem sugerida pelo 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)

Eu recebo:

     Mean        SD 
384.2283  97.0445 

e:

2. Quantiles for each variable:

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

(silencioso, longe dos valores desejados)

Para o 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)])))

um recebe:

 1000.149 0.8827613

(muito próximo dos valores desejados)

O segundo resultado está muito mais próximo dos valores reais. Mas fica pior. Se classificarmos como outliers aquelas observações para as quais o escore estimado é maior que (lembre-se de que o anterior é que F é gaussiano), a abordagem bayesiana descobrirá que todas as observações são outliers (o procedimento robusto, ao contrário, sinaliza tudo e somente os outliers como tais). Isso também implica que, se você executar uma análise bayesiana usual (não robusta) sobre os dados não classificados como outliers pelo procedimento robusto, você deve se sair bem (por exemplo, cumprir os objetivos estabelecidos em sua pergunta).$z$th$F$
Este é apenas um exemplo, mas na verdade é bastante direto mostrar que (e isso pode ser feito formalmente, veja, por exemplo, no capítulo 2 de [1]), os parâmetros de uma distribuição estudante ajustada a dados contaminados não podem depender da revelação os outliers. $t$

• [1] Ricardo A. Maronna, Douglas R. Martin, Victor J. Yohai (2006). Estatística Robusta: Teoria e Métodos (Série Wiley em Probabilidade e Estatística).
• Huber, PJ (1981). Estatísticas robustas. Nova York: John Wiley and Sons.

Na análise bayesiana, usar a distribuição gama inversa como prioridade para a precisão (o inverso da variância) é uma escolha comum. Ou a distribuição inversa do Wishart para modelos multivariados. A adição de uma prévia à variação melhora a robustez contra valores discrepantes.

Há um belo artigo de Andrew Gelman: "Distribuições anteriores para parâmetros de variação em modelos hierárquicos", onde ele discute quais boas escolhas para os anteriores sobre as variações podem ser.

$\mu$$N$$\sigma^2$$\mu$$t$$N$

$\sigma$$D$

Eu segui a discussão da pergunta original. Rasmus, quando você diz robustez, tenho certeza de que você quer dizer nos dados (discrepantes, não falta de especificação de distribuições). Tomarei a distribuição dos dados como distribuição de Laplace em vez de uma distribuição t; então, como na regressão normal em que modelamos a média, aqui modelaremos a regressão mediana (muito robusta) ou mediana (todos sabemos). Seja o modelo:

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

Obviamente, nosso objetivo é estimar os parâmetros do modelo. Esperamos que nossos anteriores sejam vagos e tenham um modelo objetivo. O modelo em mãos tem uma parte posterior da forma$f(\beta,\sigma,Y,X)$. Dando$\beta$um anterior normal com grande variância torna esse anterior vago e um anterior ao quadrado com pequenos graus de liberdade para imitar o anterior de jeffrey (anterior vago) é dado a$\sigma^2$. Com um amostrador Gibbs, o que acontece? normal anterior + aparência local = ???? nós sabemos. Também qui-quadrado anterior + probabilidade de laplace = ??? nós não sabemos a distribuição. Felizmente para nós, existe um teorema em (Aslan, 2010) que transforma a probabilidade de um lugar em uma mistura escalável de distribuições normais que nos permite desfrutar das propriedades conjugadas de nossos anteriores. Eu acho que todo o processo descrito é totalmente robusto em termos de valores discrepantes. Em um cenário multivariado, o qui-quadrado se torna uma distribuição wishart e usamos laplace multivariado e distribuições normais.

Suponha que você tenha $K$ grupos e você deseja modelar a distribuição de suas variações de amostra, talvez em relação a algumas covariáveis $\bf{x}$. Ou seja, suponha que seus dados apontem para o grupo$k \in {1 \ldots K}$ é $\textrm{Var}(y_k) \in [0, \infty)$. A pergunta aqui é: "Qual é um modelo robusto para a probabilidade de variação da amostra?" Uma maneira de abordar isso é modelar os dados transformados$\textrm{ln}[\textrm{Var}(y_k)]$ 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 \rightarrow \infty$, 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.

Um raciocínio semelhante se aplica se você estiver atribuindo uma distribuição anterior sobre um parâmetro de escala para uma distribuição normal. Tangencialmente, as distribuições lognormal e gama inversa não são aconselháveis ​​se você deseja formar um limite evitando anteriormente para fins de aproximação do modo posterior, porque eles atingem um pico acentuado se você os parametrizar para que o modo fique próximo de zero. Veja o capítulo 13 da BDA3 para discussão. Portanto, além de identificar um modelo robusto em termos de espessura da cauda, ​​lembre-se de que a curtose também é importante para sua inferência.

Espero que isso ajude você tanto quanto sua resposta a uma das minhas perguntas recentes me ajudou.