Dopasowanie rozkładu t w R: parametr skalowania

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ć?

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.

$\nu$$t$

$\nu$$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")

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

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