O que seria um modelo Bayesiano robusto para estimar a escala de uma distribuição aproximadamente normal?

32

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 σ=MAD1.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:yi

yit(m,s,ν)

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ãomsνm,sνmyisyisννNorm(μ=0,σ=1)sseria 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 
Rasmus Bååth
fonte
Talvez não seja suficientemente robusto, mas a distribuição qui-quadrado é o conjugado geralmente escolhido antes do inverso da variância.
precisa saber é o seguinte
Convém verificar se a primeira resposta a esta pergunta stats.stackexchange.com/questions/6493/… é suficiente para você; pode muito bem não ser, mas talvez seja.
jbowman
O que você tem antes do nível de contaminação? A contaminação será sistemática? Aleatória? Será gerado por uma única distribuição ou múltiplas distribuições? Sabemos algo sobre a (s) distribuição (ões) de ruído? Se pelo menos algumas das coisas acima forem conhecidas, poderíamos ajustar algum tipo de modelo de mistura. Caso contrário, não tenho certeza de quais são suas crenças sobre esse problema, e se você não tiver nenhuma, isso parecerá um cenário muito vago. Você precisa consertar algo, caso contrário, você pode escolher um ponto aleatoriamente e declarar que ele é o único ponto gerado gaussiano.
significa significado
Mas, em geral, você pode ajustar uma distribuição t que seja mais resistente a valores discrepantes ou uma mistura de distribuições t. Tenho certeza de que existem muitos trabalhos, aqui está um do Bishop research.microsoft.com/en-us/um/people/cmbishop/downloads/… e aqui está um pacote R para combinar misturas: maths.uq.edu. au / ~ gjm / mix_soft / EMMIX_R / EMMIX-manual.pdf
significa o significado
1
O seu é verdadeiro para uma população normalmente distribuída, mas não para a maioria das outras distribuiçõesσ=MAD1.4826
Henry

Respostas:

10

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:νσ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 exemplo, em , obtenho f ( ν ) = 1,18 . O estimador desejado é, em seguida, σ = s / f ( ν ) .ν=4f(ν)=1.18σ^=s/f(ν)

Tom Minka
fonte
1
Boa resposta (+1). 'no sentido de que uma única observação não pode tornar as estimativas arbitrariamente grandes', então o ponto de ruptura é 2 / n (eu estava pensando nisso) ... Como um ponto de comparação, para o procedimento ilustrado em minha resposta, é n / 2.
user603
Uau, obrigada! Pergunta de acompanhamento confusa. Será que, na verdade, faria sentido "corrigir" a escala para que ela seja consistente com o SD no caso Normal? O caso de uso em que estou pensando é ao relatar uma medida de spread. Eu não teria nenhum problema com a escala de relatórios, mas seria bom relatar algo que seria consistente com o DS, pois é a medida de disseminação mais comum (pelo menos em psicologia). Você vê uma situação em que essa correção levaria a estimativas estranhas e inconsistentes?
Rasmus Bååth
6

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νmyi

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:
    zi=|ximed(x)|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.zi>qα(z|xF)αzxFF
  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/22

a parte central dos dados é :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).zthF
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.
user603
fonte
1
Bem, o t é frequentemente proposto como uma alternativa robusta à distribuição normal. Não sei se isso é no sentido fraco ou não. Veja, por exemplo: Lange, KL, Little, RJ e Taylor, JM (1989). Modelagem estatística robusta usando a distribuição t.Jornal da Associação Estatística Americana , 84 (408), 881-896. pdf
Rasmus
1
Este é o sentido fraco. Se você possui um código R que implementa o procedimento sugerido, ficarei feliz em ilustrar minha resposta com um exemplo. caso contrário, você pode obter mais explicações no capítulo 2 deste livro.
user603
O procedimento que sugiro é basicamente descrito aqui: indiana.edu/~kruschke/BEST incluindo o código R. Vou ter que pensar na sua solução! No entanto, não parece bayesiano no sentido de que não modela todos os dados, apenas o subconjunto que "sobrevive" à etapa 2.
Rasmus Baath
1
Eu já fiz isso!
Rasmus Bååth
1

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.

jpmuc
fonte
4
Sinto muito, mas não vejo como isso responde à pergunta. Não pedi um modelo robusto antes, mas um modelo robusto .
Rasmus Bååth
0

μNσ2μtN

σD

D|μ,σN(μ,σ2)
D(d1,,dN)
p(D|μ,σ2)=1(2πσ)Nexp(N2σ2((mμ2)+s2))
ms2
m=1Ni=1Ndis2=1Ni=1Ndi2m2
p(μ,σ2|D)p(D|μ,σ2)p(μ,σ2)
(μ,σ2)p(μ,σ2|D)p(σ2|D)
σ2|DIG(α+N/2,2β+Ns2)α,β>0
σ2αβtμ
yannick
fonte
1
σ2
1
Tudo depende do que você quer dizer com robusto. O que você está dizendo agora é que deseja robustez de dados. O que eu estava propondo era robustez na especificação incorreta do modelo. Ambos são tipos diferentes de robustez.
yannick
2
Eu diria que os exemplos que dei, MAD e usando na distribuição como distribuição dos dados são exemplos de robustez em relação aos dados.
Rasmus Bååth
Eu diria que Rasmus está certo e Gelman também no BDA3, assim como um entendimento básico de que essa distribuição tem caudas mais gordas do que o normal para o mesmo parâmetro de localização.
Equilíbrio Brash
0

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=βX+ϵ, ϵ tem laplace(0 0,σ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 formaf(β,σ,Y,X). Dandoβ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σ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.

Chamberlain Foncha
fonte
2
Sua solução parece estar focada em uma estimativa robusta da localização (média / mediana). Minha pergunta foi sobre estimativa de escala com a propriedade de consistência em relação à recuperação do SD quando a distribuição de geração de dados é realmente normal.
Rasmus Bååth
Com uma estimativa robusta da localização, a escala em função da localização se beneficia imediatamente da robustez da localização. Não há outra maneira de tornar a balança robusta.
Chamberlain Foncha
De qualquer forma, devo dizer que estou esperando ansiosamente para ver como esse problema será resolvido, principalmente com uma distribuição normal, como você enfatizou.
Chamberlain Foncha 14/02
0

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 x. Ou seja, suponha que seus dados apontem para o grupok1...K é Var(yk)[0 0,). A pergunta aqui é: "Qual é um modelo robusto para a probabilidade de variação da amostra?" Uma maneira de abordar isso é modelar os dados transformadosln[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.

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.

Equilíbrio Brash
fonte
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