32

Il existe un certain nombre d' estimateurs d'échelle robustes . Un exemple notable est l’écart absolu médian qui se rapporte à l’écart type sous la forme . Dans un cadre bayésien, il existe un certain nombre de moyens pour estimer de manière fiable l' emplacement d'une distribution à peu près normale (disons une normale contaminée par des valeurs aberrantes), par exemple, on pourrait supposer que les données sont distribuées de la même manière que la distribution de Laplace. Maintenant ma question: $\sigma = \mathrm{MAD}\cdot1.4826$

Que serait un modèle bayésien permettant de mesurer de manière robuste l’ échelle d’une distribution à peu près normale, robuste au même sens que le MAD ou des estimateurs robustes similaires?

Comme dans le cas de MAD, il serait judicieux que le modèle bayésien puisse approcher le SD d'une distribution normale dans le cas où la distribution des données est réellement distribuée.

éditer 1:

Un exemple typique de modèle qui résiste à la contamination / aux valeurs aberrantes en supposant que les données est à peu près normale utilise la distribution à la manière suivante: $y_i$

y_{i} \sim t (m, s, ν)

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

Où est la moyenne, l’échelle et le degré de liberté. Avec des a priori convenables sur et , sera une estimation de la moyenne de qui sera robuste contre les valeurs aberrantes. Cependant, ne sera pas une estimation cohérente du SD de car dépend de . Par exemple, si est fixé à 4.0 et que le modèle ci-dessus est ajusté à un très grand nombre d'échantillons d'une distribution alors $m$ $s$ $\nu$ $m, s$ $\nu$ $m$ $y_i$ $s$ $y_i$ $s$ $\nu$ $\nu$ $\mathrm{Norm}(\mu=0,\sigma=1)$ $s$ serait autour de 0,82. Ce que je recherche, c’est un modèle robuste, comme le modèle t, mais pour le SD au lieu de (ou en plus) de la moyenne.

éditer 2:

Voici ci-dessous un exemple codé dans R et JAGS de la manière dont le modèle t mentionné ci-dessus est plus robuste par rapport à la moyenne.

# 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

r bayesian estimation standard-deviation robust Rasmus Bååth
la source

Peut-être que ce n'est pas assez robuste, mais la distribution chi-carré est le conjugué généralement choisi avant l'inverse de la variance.

Mike Dunlavey

Vous voudrez peut-être voir si la première réponse à cette question stats.stackexchange.com/questions/6493/… vous suffit; ce n'est peut-être pas le cas, mais c'est peut-être le cas.

jbowman

Quelle est votre priorité pour le niveau de contamination? La contamination sera-t-elle systématique? Au hasard? Sera-t-il généré par une seule distribution ou par plusieurs distributions? Savons-nous quelque chose à propos de la distribution du bruit? Si au moins certaines des choses ci-dessus sont connues, nous pourrions alors utiliser une sorte de modèle de mélange. Sinon, je ne suis pas sûr de vos croyances sur ce problème, et si vous n'en avez pas, cela vous semblera très vague. Vous devez réparer quelque chose, sinon vous pouvez choisir un point au hasard et le déclarer comme étant le seul point généré par Gauss.

signifie-à-sens

Mais en général, vous pouvez choisir une distribution t qui est plus résistante aux valeurs aberrantes ou un mélange de t-distributions. Je suis sûr qu'il existe de nombreux articles. En voici un de Bishop research.microsoft.com/en-us/um/people/cmbishop/downloads/… et voici un package R qui convient aux mélanges: maths.uq.edu. au / ~ gjm / mix_soft / EMMIX_R / EMMIX-manual.pdf

sens à la signification

1

Votre est vrai pour une population normalement distribuée, mais pas pour la plupart des autres distributions

σ = M A D \cdot 1.4826

$\sigma = \mathrm{MAD}\cdot1.4826$

Henry,

10

L'inférence bayésienne dans un modèle de bruit T avec un préalable approprié donnera une estimation robuste de l'emplacement et de l'échelle. Les conditions précises que la probabilité et le besoin préalable à satisfaire sont données dans l'étude Modélisation bayésienne de la robustesse des paramètres de localisation et d'échelle par Andrade et O'Hagan (2011). Les estimations sont robustes en ce sens qu’une seule observation ne peut pas donner des estimations arbitrairement grandes, comme le montre la figure 2 du document.

Lorsque les données sont normalement distribuées, le SD de la distribution T ajustée (pour fixe ) ne correspond pas au SD de la distribution génératrice. Mais c'est facile à résoudre. Soit l'écart-type de la distribution génératrice et soit l'écart-type de la distribution T ajustée. Si les données sont mises à l'échelle par 2, alors, sous la forme de vraisemblance, nous savons que doit être mis à l'échelle par 2. Cela implique que pour une fonction fixe . Cette fonction peut être calculée numériquement par simulation à partir d'une normale standard. Voici le code pour faire ceci: $\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)))

Par exemple, à je reçois . L'estimateur souhaité est alors . $\nu=4$ $f(\nu)=1.18$ $\hat{\sigma} = s/f(\nu)$

Tom Minka
la source

1

Bonne réponse (+1). "dans le sens où une seule observation ne peut pas rendre les estimations arbitrairement grandes", donc le point de ventilation est 2 / n (je m'interrogeais à ce sujet) ... À titre de comparaison, la procédure illustrée dans ma réponse est: n / 2.

user603

Ouah merci! Fuzzy question de suivi. Serait-il alors logique de "corriger" l’échelle pour qu’elle soit cohérente avec le DS dans le cas normal? Le cas d'utilisation auquel je songe est de signaler une mesure de propagation. Je n'aurais aucun problème avec l'échelle de rapport, mais il serait bien de signaler quelque chose qui serait compatible avec le DD car il s'agit de la mesure de propagation la plus courante (du moins en psychologie). Voyez-vous une situation où cette correction conduirait à des estimations étranges et incohérentes?

Rasmus Bååth

6

Comme vous posez une question sur un problème très précis (estimation robuste), je vais vous donner une réponse tout aussi précise. Cependant, je commencerai tout d'abord par essayer de dissiper toute hypothèse injustifiée. Il est pas vrai qu'il ya une estimation bayésienne robuste de l' emplacement (il y a des estimateurs bayésiens des endroits , mais comme je l' ai ci - dessous montrent qu'ils ne sont pas robustes et, apparemment , même le plus simple estimateur robuste de l' emplacement n'est pas bayesien). À mon avis, les raisons de l'absence de chevauchement entre les paradigmes «bayésien» et «robuste» dans le cas des localisations expliquent en grande partie pourquoi il n'existe pas non plus d'estimateurs de la dispersion à la fois robustes et bayésiens.

Avec des a priori convenables sur et , sera une estimation de la moyenne de qui sera robuste contre les valeurs aberrantes. $m, s$ $\nu$ $m$ $y_i$

En fait non. Les estimations résultantes ne seront robustes que dans un sens très faible du mot robuste. Cependant, lorsque nous disons que la médiane est robuste aux valeurs aberrantes, nous entendons le mot robuste dans un sens beaucoup plus fort. C'est-à-dire que dans les statistiques robustes, la robustesse de la médiane fait référence à la propriété que si vous calculez la médiane sur un ensemble d'observations de données provenant d'un modèle unimodal et continu, puis que vous remplacez moins de la moitié de ces observations par des valeurs arbitraires , la valeur de la médiane calculée sur les données contaminées est proche de celle que vous auriez obtenue si vous l'aviez calculée sur le jeu de données d'origine (non contaminé). Ensuite, il est facile de montrer que la stratégie d’estimation que vous proposez dans le paragraphe que j’ai cité ci-dessus n’est certainement pas robuste dans le sens où le mot est généralement compris comme la médiane.

Je ne connais pas vraiment l'analyse bayésienne. Cependant, je me demandais ce qui ne va pas dans la stratégie suivante, car elle semble simple, efficace et n’a pas encore été prise en compte dans les autres réponses. Le prior est que la bonne partie des données est tirée d'une distribution symétrique et que le taux de contamination est inférieur à la moitié. Ensuite, une stratégie simple consisterait à: $F$

calcule la médiane / folie de votre jeu de données. Puis calculez: $z_{i} = \frac{| x_{i} - med (x) |}{mad (x)}$ $z_i=\frac{|x_i-\mbox{med}(x)|}{\mbox{mad}(x)}$
exclure les observations pour lesquelles (ce qui est le quantile de la répartition des où ). Cette quantité est disponible pour de nombreux choix de et peut être initialisée pour les autres. $z_i>q_{\alpha}(z|x\sim F)$ $\alpha$ $z$ $x\sim F$ $F$
Exécuter une analyse bayésienne (habituelle, non robuste) sur les observations non rejetées.

MODIFIER:

Merci à l’opérateur d’avoir fourni un code R autonome pour effectuer une analyse bayésienne parfaite du problème.

Le code ci-dessous compare l’approche bayésienne suggérée par le PO à une alternative à la littérature de statistiques robustes (par exemple, la méthode d’ajustement proposée par Gauss pour le cas où les données peuvent contenir autant que valeurs aberrantes et la distribution des valeurs aberrantes). une bonne partie des données est gaussienne). $n/2-2$

la partie centrale des données est : $\mathcal{N}(1000,1)$

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

Ajoutez une certaine quantité de contaminants:

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

l'indice w prend la valeur 1 pour les valeurs aberrantes. Je commence par l'approche proposée par le PO:

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)

Je reçois:

     Mean        SD 
384.2283  97.0445

et:

2. Quantiles for each variable:

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

(calme loin donc des valeurs cibles)

Pour la méthode robuste,

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

on obtient:

 1000.149 0.8827613

(très proche des valeurs cibles)

$z$ th $F$
$t$

[1] Ricardo A. Maronna, Douglas R. Martin et Victor J. Yohai (2006). Statistiques robustes: théorie et méthodes (série de Wiley sur les probabilités et les statistiques).
Huber, PJ (1981). Statistiques robustes. New York: John Wiley et ses fils.

utilisateur603
la source

1

Eh bien, le t est souvent proposé comme une alternative robuste à la distribution normale. Je ne sais pas si c'est dans le sens faible ou pas. Voir, par exemple: Lange, KL, Little, RJ et Taylor, JM (1989). Modélisation statistique robuste utilisant la distribution t. Journal de l'Association américaine de statistique , 84 (408), 881-896. pdf

Rasmus Bååth

1

C'est le sens faible. Si vous avez un code R qui implémente la procédure que vous suggérez, je me ferai un plaisir d’illustrer ma réponse avec un exemple. sinon, vous pouvez obtenir plus d'explications au chapitre 2 de ce manuel.

user603

La procédure que je suggère est essentiellement décrite ici: indiana.edu/~kruschke/BEST, y compris le code R. Je vais devoir réfléchir à votre solution! Il ne semble toutefois pas être bayésien dans le sens où il ne modélise pas toutes les données, mais seulement le sous-ensemble qui "survit" à l'étape 2.

Rasmus Bååth

\to \infty

$\rightarrow \infty$

1

Je l'ai maintenant fait!

Rasmus Bååth

1

Dans l'analyse bayésienne, l'utilisation de la distribution gamma inverse en tant qu'antérieur de la précision (l'inverse de la variance) est un choix courant. Ou la distribution Wishart inverse pour les modèles multivariés. L'ajout d'un prior sur la variance améliore la robustesse contre les valeurs aberrantes.

Andrew Gelman a rédigé un article intéressant intitulé "Distributions préalables des paramètres de variance dans les modèles hiérarchiques", dans lequel il explique quels sont les bons choix pour les variables a priori sur les variances.

jpmuc
la source

4

Je suis désolé mais je ne vois pas comment cela répond à la question. Je n'ai pas demandé un modèle robuste, mais plutôt un modèle robuste .

Rasmus Bååth le

0

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

$\sigma$ $D$

{D |}_{μ, σ} \sim N (μ, σ^{2})

$\left.D\right|_{\mu,\sigma} \sim \mathcal{N}(\mu,\sigma^2)$

D \equiv (d_{1}, \dots, d_{N})

$D \equiv (d_1,\ldots,d_N)$

p (D | μ, σ^{2}) = \frac{1}{(\sqrt{2 π} σ)^{N}} \exp (- \frac{N}{2 σ^{2}} ((m - μ^{2}) + s^{2}))

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

m

$m$

s^{2}

$s^2$

m = \frac{1}{N} \sum_{i = 1}^{N} d_{i} s^{2} = \frac{1}{N} \sum_{i = 1}^{N} d_{i}^{2} - m^{2}

$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$

p (μ, σ^{2} | D) \propto p (D | μ, σ^{2}) p (μ, σ^{2})

$p(\mu,\sigma^2|D) \propto p(D|\mu,\sigma^2) p(\mu,\sigma^2)$

(μ, σ^{2})

$(\mu,\sigma^2)$

p (μ, σ^{2} | D)

$p(\mu,\sigma^2|D)$

p (σ^{2} | D)

$p(\sigma^2|D)$

{σ^{2} |}_{D} \sim I G (α + N / 2, 2 β + N s^{2}) α, β > 0

$\left.\sigma^2\right|_{D} \sim \mathcal{IG}\left(\alpha+N/2,2\beta+Ns^2\right) \qquad \alpha,\beta>0$

σ^{2}

$\sigma^2$

α

$\alpha$

β

$\beta$

t

$t$

μ

$\mu$

yannick
la source

1

σ^{2}

$σ^2$

1

Tout dépend de ce que vous entendez par robuste. Ce que vous dites en ce moment, c'est que vous aimeriez la robustesse par rapport aux données. Ce que je proposais était la robustesse par rapport à une spécification erronée du modèle. Ce sont deux types différents de robustesse.

Yannick

2

Je dirais que les exemples que j'ai donnés, MAD et utilisant at distribution comme distribution des données sont des exemples de robustesse vis-à-vis des données.

Rasmus Bååth

Je dirais que Rasmus a raison, tout comme Gelman dans BDA3, de même qu'une compréhension de base selon laquelle cette distribution a des queues plus épaisses que la normale pour le même paramètre d'emplacement

Brash Equilibrium

0

J'ai suivi la discussion de la question initiale. Rasmus quand vous parlez de robustesse, je suis sûr que vous voulez dire dans les données (valeurs aberrantes, pas une spécification erronée des distributions). Je prendrai la distribution des données comme étant la distribution de Laplace au lieu d’une distribution t, puis, comme dans la régression normale dans laquelle nous modélisons la moyenne, nous modélisons ici la régression médiane (très robuste), également connue (connue de tous). Que le modèle soit:

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

$f(\beta,\sigma,Y,X)$ $\beta$ $\sigma^2$ . Que se passe-t-il avec un échantillonneur Gibbs? avant normal + laplace likehood = ???? nous savons. Aussi, khi-deux avant + laplace vraisemblance = ??? nous ne connaissons pas la distribution. Heureusement pour nous, il existe un théorème dans (Aslan, 2010) qui transforme une vraisemblance en un mélange à l'échelle de distributions normales qui nous permet ensuite de profiter des propriétés conjuguées de nos a priori. Je pense que l'ensemble du processus décrit est totalement robuste en termes de valeurs aberrantes. Dans un paramètre multivarié, le chi-carré devient une distribution wishart et nous utilisons des distributions multivariées laplace et normales.

Chamberlain Foncha
la source

2

Votre solution semble être axée sur une estimation robuste de la position (moyenne / médiane). Ma question concernait plutôt l’estimation d’échelle avec la propriété de cohérence vis-à-vis de la récupération de la SD lorsque la distribution génératrice de données est en réalité normale.

Rasmus Bååth

Avec une estimation robuste de la localisation, l’échelle en fonction de la localisation bénéficie immédiatement de la robustesse de la localisation. Il n'y a pas d'autre moyen de rendre la balance robuste.

Chamberlain Foncha

Quoi qu'il en soit, je dois dire que j'attends avec impatience de voir comment ce problème sera résolu plus particulièrement avec une distribution normale, comme vous l'avez souligné.

Chamberlain Foncha

0

$K$ $\bf{x}$ $k \in {1 \ldots K}$ $\textrm{Var}(y_k) \in [0, \infty)$ $\textrm{ln}[\textrm{Var}(y_k)]$ $t$ $n \rightarrow \infty$

Un raisonnement similaire s'applique si, au lieu de cela, vous affectez une distribution antérieure à un paramètre d'échelle pour une distribution normale. Tangentiellement, les distributions log-normales et gamma-inverses ne sont pas recommandées si vous souhaitez former une limite en évitant le précédent aux fins de l'approximation en mode postérieur, car elles atteignent des sommets nets si vous les paramétrez de sorte que le mode soit proche de zéro. Voir le chapitre 13 de BDA3 pour la discussion. Donc, en plus d'identifier un modèle robuste en termes d'épaisseur de la queue, gardez à l'esprit que le kurtosis peut également être important pour votre inférence.

J'espère que cela vous aidera autant que votre réponse à l'une de mes dernières questions m'a aidé.

Équilibre de Brash
la source

1

Ma question portait sur la situation lorsque vous avez un groupe et sur la façon d’évaluer solidement l’échelle de ce groupe. Dans le cas des valeurs aberrantes, je ne crois pas que la variance de l'échantillon soit considérée comme robuste.

Rasmus Bååth

Si vous avez un groupe et que vous estimez sa distribution normale, votre question s'applique alors à la forme du paramètre précédent sur son paramètre d'échelle. Comme ma réponse le laisse supposer, vous pouvez utiliser la distribution pour la transformation de son journal ou choisir une distribution à la queue épaisse avec un support réel positif, en faisant attention aux autres aspects de cette distribution, tels que le kurtosis. En bout de ligne, si vous recherchez un modèle robuste pour un paramètre d'échelle, utilisez-le à la distribution par rapport à sa transformation de journal ou à une autre distribution grasse.

Brash Equilibrium

Que serait un modèle bayésien robuste pour estimer l’échelle d’une distribution approximativement normale?

Réponses:

MODIFIER: