Istnieje wiele niezawodnych estymatorów skali . Godnym uwagi przykładem jest mediana bezwzględnego odchylenia, które odnosi się do odchylenia standardowego jako $\sigma = \mathrm{MAD}\cdot1.4826$ . W ramach bayesowskich istnieje wiele sposobów dokładnego oszacowania lokalizacji mniej więcej normalnej dystrybucji (powiedzmy normalnej zanieczyszczonej wartościami odstającymi), na przykład można założyć, że dane są dystrybuowane tak jak w dystrybucji lub dystrybucji Laplace'a. Teraz moje pytanie:

Czym byłby model bayesowski do pomiaru skali mniej więcej normalnego rozkładu w solidny sposób, solidny w tym samym sensie, co MAD lub podobne niezawodne estymatory?

Podobnie jak w przypadku MAD, byłoby fajnie, gdyby model bayesowski mógł zbliżyć się do SD rozkładu normalnego w przypadku, gdy rozkład danych jest faktycznie rozkładem normalnym.

edycja 1:

Typowym przykładem modelu, który jest odporny na zanieczyszczenia / skrajnych przy założeniu danych $y_i$ jest w przybliżeniu normalnie stosuje się rozkład, takich jak:

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

Gdzie $m$ jest średnią, $s$ jest skalą, a $\nu$ jest stopniem swobody. Z odpowiednich priors na $m, s$ i $\nu$ , $m$ będą szacunkową średnią $y_i$ który będzie odporny na błędne. Jednak $s$ nie będzie spójne oszacowanie SD $y_i$ jak $s$ zależy $\nu$ . Na przykład, jeśli $\nu$ będzie ustalona na 4,0 i wzór powyżej, być przymocowane do ogromnej ilości próbek z rozkład to s będzie wynosić około 0,82. To, czego szukam, to model, który jest solidny, podobnie jak model t, ale dla SD zamiast (lub dodatkowo) średniej.$\mathrm{Norm}(\mu=0,\sigma=1)$$s$

edycja 2:

Poniżej znajduje się zakodowany przykład w R i JAGS, w którym wspomniany powyżej model t jest bardziej niezawodny w odniesieniu do średniej.

# 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 

Wnioskowanie bayesowskie w modelu hałasu T z odpowiednim wyprzedzeniem da wiarygodne oszacowanie lokalizacji i skali. Dokładne warunki, które prawdopodobieństwo i wcześniejsze potrzeby muszą spełnić, podane są w pracy Bayesowskiego modelowania odporności parametrów lokalizacji i skali przez Andrade'a i O'Hagana (2011). Szacunki są wiarygodne w tym sensie, że pojedyncza obserwacja nie może uczynić oszacowań arbitralnie dużymi, jak pokazano na ryc. 2 artykułu.

Kiedy dane są normalnie dystrybuowane, SD dopasowanego rozkładu T (dla stałego ) nie pasuje do SD rozkładu generującego. Ale łatwo to naprawić. Niech σ być odchylenie standardowe rozkładu prądotwórczego i niech s jest odchyleniem standardowym rozkładu dopasowanego T. Jeżeli dane są skalowane o 2, to z postaci prawdopodobieństwa wiemy, że s musi być skalowane o 2. Oznacza to, że s = σ f ( ν ) dla niektórych stałych funkcji f . Ta funkcja może być obliczona numerycznie przez symulację standardowej normy. Oto kod, aby to zrobić:$\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)))

Na przykład przy otrzymuję f ( ν ) = 1,18 . Żądany estymator następnie σ = s / f ( ν ) .$\nu=4$$f(\nu)=1.18$$\hat{\sigma} = s/f(\nu)$

Kiedy zadajesz pytanie na temat bardzo precyzyjnego problemu (rzetelne oszacowanie), dam ci równie precyzyjną odpowiedź. Najpierw jednak zacznę próbować rozwiać nieuzasadnione założenie. Nie jest prawdą, że istnieje solidna bayesowska ocena lokalizacji (istnieją bayesowskie estymatory lokalizacji, ale jak ilustruję poniżej, nie są one solidne i, wydaje , nawet najprostszy solidny estymator lokalizacji nie jest bayesowski). Moim zdaniem powody braku nakładania się paradygmatu „bayesowskiego” i „solidnego” w przypadku lokalizacji znacznie przyczyniają się do wyjaśnienia, dlaczego nie ma również estymatorów rozproszenia, które byłyby zarówno solidne, jak i bayesowskie.

Z odpowiednich priors na i v , m będą szacunkową średnią y $m, s$$\nu$$m$$y_i$ który będzie odporny na błędne.

Właściwie nie. Wynikowe oszacowania będą solidne tylko w bardzo słabym znaczeniu tego słowa. Jednak gdy mówimy, że mediana jest solidna na wartości odstające, mamy na myśli słowo solidne w znacznie silniejszym znaczeniu. To znaczy, w solidnych statystykach, odporność mediany odnosi się do właściwości, która jeśli obliczymy medianę na zbiorze danych obserwacji pochodzących z jednomodalnego modelu ciągłego, a następnie zastąpimy mniej niż połowę tych obserwacji arbitralnymi wartościami , wartość mediany obliczonej na zanieczyszczonych danych jest zbliżona do wartości, którą miałbyś, gdybyś obliczył ją na oryginalnym (niezanieczyszczonym) zbiorze danych. Łatwo zatem wykazać, że strategia szacowania zaproponowana w cytowanym powyżej akapicie zdecydowanie nie jest solidny w tym sensie, w jaki to słowo jest zazwyczaj rozumiane jako mediana.

Zupełnie nie znam analizy bayesowskiej. Zastanawiałem się jednak, co jest nie tak z następującą strategią, ponieważ wydaje się ona prosta, skuteczna, a jednak nie została uwzględniona w innych odpowiedziach. Pierwszym z nich jest to, że duża część danych pochodzi z symetrycznego rozkładu i że stopień zanieczyszczenia jest mniejszy niż połowa. Następnie prostą strategią byłoby:$F$

1. obliczyć medianę / szalenie z zestawu danych. Następnie oblicz: $$z_i=\frac{|x_i-\mbox{med}(x)|}{\mbox{mad}(x)}$$
2. wykluczyć obserwacje, dla których (jest to kwantyl α rozkładu z, gdy x ∼ F ). Ta ilość jest dostępna dla wielu opcji wyboru $z_i>q_{\alpha}(z|x\sim F)$$\alpha$$z$$x\sim F$$F$ i może zostać załadowana do innych.
3. Przeprowadź (zwykłą, mało wiarygodną) analizę bayesowską na nie odrzuconych obserwacjach.

EDYTOWAĆ:

Dzięki OP za udostępnienie samodzielnego kodu R do przeprowadzenia bayesowskiej analizy problemu bonna fide.

poniższy kod porównuje podejście bayesowskie sugerowane przez PO z jego alternatywą dla solidnej literatury statystycznej (np. metoda dopasowania zaproponowana przez Gaussa w przypadku, gdy dane mogą zawierać nawet $n/2-2$ skrajnych i dystrybucji duża część danych to Gaussa).

centralna część danych to :$\mathcal{N}(1000,1)$

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

Dodaj pewną ilość zanieczyszczeń:

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

indeks w przyjmuje wartość 1 dla wartości odstających. Zaczynam od podejścia zaproponowanego przez 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)

Dostaję:

     Mean        SD 
384.2283  97.0445 

i:

2. Quantiles for each variable:

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

(cicho daleko od wartości docelowych)

W przypadku niezawodnej metody

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

dostaje się:

 1000.149 0.8827613

(bardzo blisko wartości docelowych)

Drugi wynik jest znacznie bliższy rzeczywistym wartościom. Ale robi się coraz gorzej. Jeśli będziemy klasyfikować jako odstających te obserwacje, których szacunkowa -score jest większy niż (należy pamiętać, że przed to, że F jest Gaussa) wtedy Bayesa znaleziska podejście, że$z$th$F$ wszystkie obserwacje są odstające (solidna procedura, w przeciwieństwie do tego, flagi i wszystko tylko wartości odstające jako takie). Oznacza to również, że jeśli miałbyś przeprowadzić zwykłą (niesolidną) analizę bayesowską na danych niesklasyfikowanych jako odstające według solidnej procedury, powinieneś zrobić dobrze (np. Osiągnąć cele określone w pytaniu).
To tylko przykład, ale to rzeczywiście dość proste, aby pokazać, że (a może to zrobić formalnie, patrz na przykład, w rozdziale 2 [1]) parametry studenta dystrybucji dopasowane do danych skażonych nie można zależało, aby odsłonić wartości odstające. $t$

• [1] Ricardo A. Maronna, Douglas R. Martin, Victor J. Yohai (2006). Solidne statystyki: teoria i metody (seria Wileya w prawdopodobieństwie i statystyce).
• Huber, PJ (1981). Solidne statystyki. Nowy Jork: John Wiley and Sons.

W analizie bayesowskiej powszechnym wyborem jest odwrotny rozkład gamma jako pierwszeństwo dla precyzji (odwrotność wariancji). Lub odwrotny rozkład Wishart dla modeli wielowymiarowych. Dodanie wcześniejszego wariantu poprawia odporność na wartości odstające.

Jest ładny artykuł Andrew Gelmana: „Wcześniejsze rozkłady parametrów wariancji w modelach hierarchicznych”, w których omawia, jakie mogą być dobre wybory dla priorów w wariancjach.

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

$\sigma$$D$

$$\left.D\right|_{\mu,\sigma} \sim \mathcal{N}(\mu,\sigma^2)$$ $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)$$ $m$$s^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(\mu,\sigma^2|D) \propto p(D|\mu,\sigma^2) p(\mu,\sigma^2)$$ $(\mu,\sigma^2)$$p(\mu,\sigma^2|D)$$p(\sigma^2|D)$ $$\left.\sigma^2\right|_{D} \sim \mathcal{IG}\left(\alpha+N/2,2\beta+Ns^2\right) \qquad \alpha,\beta>0$$ $\sigma^2$. This estimator will be more or less tolerant to small excursions from misspecifications on the model by varying $\alpha$ and/or $\beta$. 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 $\mu$ that you mention.

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=\beta X+\epsilon$, $\epsilon$ has laplace$(0,$$\sigma^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(\beta,\sigma,Y,X)$. Giving $\beta$ 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 $\sigma^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.

Suppose that you have $K$ groups and you want to model the distribution of their sample variances, perhaps in relation to some covariates $\bf{x}$. That is, suppose that your data point for group $k \in {1 \ldots K}$ is $\textrm{Var}(y_k) \in [0, \infty)$. 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 $\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.

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.