Dopasowanie rozkładu t w R: parametr skalowania

17

Jak dopasować parametry rozkładu t, tj. Parametry odpowiadające „średniej” i „odchyleniu standardowemu” rozkładu normalnego. Zakładam, że są one nazywane „średnimi” i „skalowaniem / stopniami swobody” dla rozkładu t?

Poniższy kod często powoduje błędy „nieudana optymalizacja”.

library(MASS)
fitdistr(x, "t")

Czy najpierw muszę skalować x, czy przeliczać na prawdopodobieństwa? Jak najlepiej to zrobić?

użytkownik12719
źródło
2
Nie udaje się nie dlatego, że trzeba skalować parametry, ale dlatego, że nie działa optymalizator. Zobacz moją odpowiedź poniżej.
Sergey Bushmanov

Odpowiedzi:

16

fitdistrwykorzystuje techniki największego prawdopodobieństwa i optymalizacji w celu znalezienia parametrów danej dystrybucji. Czasami, szczególnie dla t-dystrybucji, jak zauważył @ user12719, optymalizacja w postaci:

fitdistr(x, "t")

kończy się błędem.

W takim przypadku powinieneś pomóc optymalizatorowi, podając punkt początkowy i dolną granicę, aby rozpocząć wyszukiwanie optymalnych parametrów:

fitdistr(x, "t", start = list(m=mean(x),s=sd(x), df=3), lower=c(-1, 0.001,1))

Uwaga: df=3najlepiej zgadnąć, co dfmoże być „optymalne” . Po podaniu tych dodatkowych informacji Twój błąd zniknie.

Kilka fragmentów, które pomogą Ci lepiej zrozumieć wewnętrzną mechanikę fitdistr:

W przypadku rozkładu normalnego, log-normalnego, geometrycznego, wykładniczego i Poissona stosowane są zamknięte MLE (i dokładne błędy standardowe) i startnie należy ich podawać.

...

Dla następujących nazwanych rozkładów, rozsądne wartości początkowe zostaną obliczone, jeśli startzostaną pominięte lub tylko częściowo określone: ​​„cauchy”, „gamma”, „logistic”, „ujemny dwumianowy” (sparametryzowany przez mu i rozmiar), „t” i „weibull” „. Należy zauważyć, że te wartości początkowe mogą nie być wystarczająco dobre, jeśli dopasowanie jest słabe: w szczególności nie są one odporne na wartości odstające, chyba że dopasowany rozkład jest długi.

Siergiej Bushmanow
źródło
1
Obie odpowiedzi (Flom i Bushmanov) są pomocne. Wybieram ten, ponieważ czyni to bardziej wyraźnym, że przy odpowiednich wartościach początkowych i ograniczeniach optymalizacja „fitdistr” jest zbieżna.
user12719,
10

νt

νt

set.seed(1234)
n <- 10
x <- rt(n,  df=2.5)

make_loglik  <-  function(x)
    Vectorize( function(nu) sum(dt(x, df=nu,  log=TRUE)) )

loglik  <-  make_loglik(x)
plot(loglik,  from=1,  to=100,  main="loglikelihood function for df     parameter", xlab="degrees of freedom")
abline(v=2.5,  col="red2")

wprowadź opis zdjęcia tutaj

n

Wypróbujmy niektóre symulacje:

t_nu_mle  <-  function(x) {
    loglik  <-  make_loglik(x)
    res  <-  optimize(loglik, interval=c(0.01, 200), maximum=TRUE)$maximum
    res   
}

nus  <-  replicate(1000, {x <- rt(10, df=2.5)
    t_nu_mle(x) }, simplify=TRUE)

> mean(nus)
[1] 45.20767
> sd(nus)
[1] 78.77813

Wyświetlanie oszacowania jest bardzo niestabilne (patrząc na histogram, znaczna część oszacowanych wartości znajduje się w górnej granicy podanej w celu optymalizacji 200).

Powtarzanie z większą próbką:

nus  <-  replicate(1000, {x <- rt(50, df=2.5)
    t_nu_mle(x) }, simplify=TRUE)
> mean(nus)
[1] 4.342724
> sd(nus)
[1] 14.40137

co jest znacznie lepsze, ale średnia wciąż znacznie przewyższa prawdziwą wartość 2,5.

Pamiętaj, że jest to uproszczona wersja prawdziwego problemu, w którym należy również oszacować parametry lokalizacji i skali.

tν

kjetil b halvorsen
źródło
5
Twój wniosek, że problemy z oszacowaniem df mogą faktycznie działać wbrew powodowi wyboru rozkładu t w pierwszej kolejności (tj. Odporności), wydaje się prowokować.
użytkownik12719,
1
(+1) „Bez ograniczeń powyżej” nie jest złą odpowiedzią i może być przydatne do niektórych celów w połączeniu z oszacowaniem przedziału. Ważne jest, aby nie ślepo wykorzystywać zaobserwowanych informacji Fishera do tworzenia przedziałów ufności Walda.
Scortchi - Przywróć Monikę
8

W pomocy dla fitdistr znajduje się ten przykład:

fitdistr(x2, "t", df = 9)

wskazując, że potrzebujesz tylko wartości df. Ale to zakłada standaryzację.

Dla większej kontroli pokazują również

mydt <- function(x, m, s, df) dt((x-m)/s, df)/s
fitdistr(x2, mydt, list(m = 0, s = 1), df = 9, lower = c(-Inf, 0))

gdzie parametrami byłoby m = średnia, s = odchylenie standardowe, df = stopnie swobody

Peter Flom - Przywróć Monikę
źródło
1
Chyba jestem zdezorientowany parametrami rozkładu T. Czy ma parametry 2 (średnia, df) czy 3 (średnia, odchylenie standardowe, df)? Zastanawiałem się, czy można dopasować parametr „df”.
user12719,
1
@ user12719 Rozkład t-Studenta ma trzy parametry: położenie, skalę i stopnie swobody. Nie są one określane jako średnia, odchylenie standardowe i df, ponieważ średnia i wariancja tego rozkładu zależą od trzech parametrów. W niektórych przypadkach nie istnieją. Peter Flom naprawia df, ale można to również uznać za nieznany parametr.
1
@PeterFlom W przypadku rozkładu Cauchy'ego wyraźnie widać, że m i s to lokalizacja i skala. Zgadzam się, że notacja mi sugeruje, że reprezentują one odpowiednio średnią i odchylenie standardowe. Ale to może być tylko uproszczenie \mui \sigmajak dobrze. Nawiasem mówiąc, dawno +1.
1
@PeterFlom Czy to cytowanie z pliku pomocy R sugeruje, że df ma zawsze 9 dla dystrybucji studentów? Czy nie sądzisz, że df również należy szacować? Właściwie to brak dfjest przyczyną błędu, a właściwa odpowiedź powinna dostarczyć przepis na jego znalezienie.
Sergey Bushmanov,
1
@PeterFlom BTW, jeśli przeczytasz plik pomocy kilka linii powyżej cytatu, dowiesz się, dlaczego df=9jest dobry w ich przykładzie i nie ma znaczenia tutaj.
Sergey Bushmanov,