Czym byłby solidny model bayesowski do szacowania skali mniej więcej normalnego rozkładu?

32

Istnieje wiele niezawodnych estymatorów skali . Godnym uwagi przykładem jest mediana bezwzględnego odchylenia, które odnosi się do odchylenia standardowego jako σ=MAD1.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 yi jest w przybliżeniu normalnie stosuje się rozkład, takich jak:

yit(m,s,ν)

Gdzie m jest średnią, s jest skalą, a ν jest stopniem swobody. Z odpowiednich priors na m,s i ν , m będą szacunkową średnią yi który będzie odporny na błędne. Jednak s nie będzie spójne oszacowanie SD yi jak s zależy ν . Na przykład, jeśli ν 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.Norm(μ=0,σ=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 
Rasmus Bååth
źródło
Może nie jest wystarczająco mocny, ale rozkład chi-kwadrat jest zwykle wybieranym koniugatem dla odwrotności wariancji.
Mike Dunlavey
Możesz sprawdzić, czy pierwsza odpowiedź na to pytanie stats.stackexchange.com/questions/6493/... jest dla Ciebie wystarczająca; może nie być, ale może tak jest.
łucznik
Jaki jest Twój priorytet dla poziomu zanieczyszczenia? Czy zanieczyszczenie będzie systematyczne? Losowy? Czy będzie generowany przez jedną dystrybucję, czy wiele dystrybucji? Czy wiemy coś o rozkładach hałasu? Jeśli przynajmniej niektóre z powyższych rzeczy są znane, moglibyśmy dopasować jakiś model mieszanki. W przeciwnym razie nie jestem pewien, jakie są twoje przekonania na temat tego problemu, a jeśli go nie masz, wydaje się to bardzo niejasnym ustawieniem. Musisz coś naprawić, w przeciwnym razie możesz losowo wybrać punkt i zadeklarować go jako jedyny punkt generowany przez Gaussa.
oznacza, co oznacza
Ale ogólnie można dopasować rozkład t, który jest bardziej odporny na wartości odstające, lub mieszankę rozkładów t. Jestem pewien, że jest wiele artykułów, tutaj jest jeden autorstwa Bishop research.microsoft.com/en-us/um/people/cmbishop/downloads/... a oto pakiet R pasujący do mieszanin: maths.uq.edu. au / ~ gjm / mix_soft / EMMIX_R / EMMIX-manual.pdf
znaczy znaczy
1
Twoje jest prawdziwe dla populacji normalnie rozmieszczonej, ale nie dla większości innych dystrybucjiσ=MAD1.4826
Henry

Odpowiedzi:

10

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ć:νσ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)))

Na przykład przy otrzymuję f ( ν ) = 1,18 . Żądany estymator następnie σ = s / f ( ν ) .ν=4f(ν)=1.18σ^=s/f(ν)

Tom Minka
źródło
1
Dobra odpowiedź (+1). „w tym sensie, że pojedyncza obserwacja nie może uczynić oszacowań arbitralnie dużymi”, więc punkt podziału wynosi 2 / n (zastanawiałem się nad tym) ... Dla porównania, dla procedury przedstawionej w mojej odpowiedzi jest to n / 2.
user603
Wow, dzięki! Rozmyte pytanie kontrolne. Czy zatem sensowne byłoby „skorygowanie” skali, aby była zgodna z SD w normalnym przypadku? Przypadek użycia, o którym myślę, to zgłoszenie miary rozprzestrzeniania się. Nie miałbym problemu ze skalą raportowania, ale byłoby miło zgłosić coś, co byłoby zgodne z SD, ponieważ jest to najczęstsza miara rozprzestrzeniania się (przynajmniej w psychologii). Czy widzisz sytuację, w której ta korekta prowadziłaby do dziwnych i niespójnych szacunków?
Rasmus Bååth
6

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 Im,sνmyi 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:
    zi=|ximed(x)|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 F.zi>qα(z|xF)αzxFF 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/22 skrajnych i dystrybucji duża część danych to Gaussa).

centralna część danych to :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, żezthF 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.
użytkownik603
źródło
1
Cóż, t jest często proponowane jako solidna alternatywa dla rozkładu normalnego. Nie wiem, czy to w słabym znaczeniu, czy nie. Patrz na przykład: Lange, KL, Little, RJ i Taylor, JM (1989). Solidne modelowanie statystyczne z wykorzystaniem rozkładu t. Journal of the American Statistics Association , 84 (408), 881-896. pdf
Rasmus Bååth
1
To jest słaby sens. Jeśli masz kod R, który implementuje sugerowaną przez Ciebie procedurę, chętnie zilustruję moją odpowiedź przykładem. w przeciwnym razie możesz uzyskać więcej wyjaśnień w rozdziale 2 tego podręcznika.
user603
Procedura, którą proponuję, została w zasadzie opisana tutaj: indiana.edu/~kruschke/BEST wraz z kodem R. Będę musiał pomyśleć o twoim rozwiązaniu! Nie wydaje się jednak Bayesowskie w tym sensie, że nie modeluje wszystkich danych, tylko podzbiór, który „przetrwa” krok 2.
Rasmus Bååth
1
Teraz to zrobiłem!
Rasmus Bååth
1

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.

jpmuc
źródło
4
Przykro mi, ale nie rozumiem, jak to odpowiada na pytanie. Nie prosiłem o solidny wcześniej, ale raczej o solidny model .
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. 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
źródło
1
„Solidny estymator parametru lokalizacji μ jakiegoś zbioru danych o rozmiarze N jest uzyskiwany, gdy ktoś przypisuje Jeffreysa przed wariancją σ2)rozkładu normalnego. ”Czy to nie jest normalny model, który opisujesz typowy przykład niestabilnego modelu? To znaczy, jedna wyłączona wartość może mieć duży wpływ na parametry modelu. Istnieje duża różnica między tylna nad średnią jest rozkładem t (jak w twoim przypadku), a rozkład dla danych jest rozkładem t (jak to jest powszechny przykład solidnego modelu bayesowskiego do szacowania średniej).
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
Powiedziałbym, że podane przeze mnie przykłady, MAD i używane przy dystrybucji jako dystrybucja danych są przykładami niezawodności w odniesieniu do danych.
Rasmus Bååth,
Powiedziałbym, że Rasmus ma rację, podobnie Gelman i inni w BDA3, podobnie jak podstawowe zrozumienie, że rozkład ten ma grubsze ogony niż normalny dla tego samego parametru lokalizacji
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
źródło
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
W każdym razie muszę powiedzieć, że z niecierpliwością czekam, aby zobaczyć, jak ten problem zostanie rozwiązany najbardziej, zwłaszcza przy normalnym rozkładzie, jak podkreślono.
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
źródło
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