Anpassung der t-Verteilung in R: Skalierungsparameter

17

Wie passe ich die Parameter einer t-Verteilung an, dh die Parameter, die dem Mittelwert und der Standardabweichung einer Normalverteilung entsprechen? Ich nehme an, sie heißen 'Mittelwert' und 'Skalierung / Freiheitsgrade' für eine t-Verteilung.

Der folgende Code führt häufig zu Fehlern bei der Optimierung.

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

Muss ich x zuerst skalieren oder in Wahrscheinlichkeiten konvertieren? Wie geht das am besten?

user12719
quelle
2
Dies schlägt fehl, nicht weil Sie Parameter skalieren müssen, sondern weil der Optimierer fehlschlägt. Siehe meine Antwort unten.
Sergey Bushmanov

Antworten:

16

fitdistrVerwendet Maximum-Likelihood- und Optimierungstechniken, um Parameter einer bestimmten Verteilung zu finden. Manchmal, besonders für T-Distribution, wie @ user12719 bemerkte, die Optimierung in der Form:

fitdistr(x, "t")

schlägt mit einem Fehler fehl.

In diesem Fall sollten Sie dem Optimierer eine Hand geben, indem Sie den Startpunkt und die Untergrenze angeben, um mit der Suche nach optimalen Parametern zu beginnen:

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

Beachten Sie, df=3ist Ihre beste Vermutung, was ein "Optimum" sein dfkönnte. Nachdem Sie diese zusätzlichen Informationen eingegeben haben, ist Ihr Fehler behoben.

Einige Auszüge zum besseren Verständnis der inneren Mechanik von fitdistr:

Für die Normal-, Log-Normal-, Geometrie-, Exponential- und Poisson-Verteilungen werden MLEs in geschlossener Form (und genaue Standardfehler) verwendet und startsollten nicht angegeben werden.

...

Für die folgenden benannten Verteilungen werden sinnvolle Startwerte berechnet, wenn sie startweggelassen oder nur teilweise angegeben werden: "cauchy", "gamma", "logistic", "negative binomial" (parametrisiert durch mu und size), "t" und "weibull" ". Beachten Sie, dass diese Ausgangswerte bei schlechter Passform möglicherweise nicht gut genug sind. Insbesondere sind sie nur dann gegen Ausreißer beständig, wenn die angepasste Verteilung einen langen Schwanz aufweist.

Sergey Bushmanov
quelle
1
Beide Antworten (Flom und Bushmanov) sind hilfreich. Ich wähle dieses aus, weil es expliziter macht, dass mit den richtigen Anfangswerten und Einschränkungen die 'fitdistr'-Optimierung konvergiert.
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")

Bildbeschreibung hier eingeben

n

Lassen Sie uns einige Simulationen versuchen:

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

Das Anzeigen der Schätzung ist sehr instabil (im Histogramm befindet sich ein beträchtlicher Teil der geschätzten Werte an der Obergrenze für die Optimierung von 200).

Wiederholen mit einer größeren Stichprobe:

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

Das ist viel besser, aber der Mittelwert liegt immer noch weit über dem wahren Wert von 2,5.

Denken Sie dann daran, dass dies eine vereinfachte Version des eigentlichen Problems ist, bei dem auch Standort- und Skalenparameter geschätzt werden müssen.

tν

kjetil b halvorsen
quelle
5
Ihre Schlussfolgerung, dass die Probleme bei der Schätzung von df tatsächlich gegen den Grund für die Wahl einer t-Verteilung (dh Robustheit) wirken könnten, ist nachdenklich.
user12719
1
(+1) "Ungebunden über" ist keine falsche Antwort und kann in Verbindung mit einer Intervallschätzung für einige Zwecke nützlich sein. Es ist wichtig, die beobachteten Fisher-Informationen nicht blind zu verwenden, um Wald-Vertrauensintervalle zu bilden.
Scortchi
8

In der Hilfe für fitdistr ist dieses Beispiel:

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

Zeigt an, dass Sie nur einen Wert für df benötigen. Das setzt aber eine Standardisierung voraus.

Für mehr Kontrolle zeigen sie auch

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

wobei die Parameter m = Mittelwert sein würden, s = Standardabweichung, df = Freiheitsgrade

Peter Flom - Wiedereinsetzung von Monica
quelle
1
Ich glaube, ich bin verwirrt über die Parameter einer t-Verteilung. Hat es 2 (Mittelwert, df) oder 3 (Mittelwert, Standardabweichung, df) Parameter? Ich habe mich gefragt, ob man den Parameter 'df' anpassen kann.
user12719
1
@ user12719 Die Student's-t-Verteilung enthält drei Parameter: Position, Maßstab und Freiheitsgrade. Sie werden nicht als Mittelwert, Standardabweichung und df bezeichnet, da der Mittelwert und die Varianz dieser Verteilung von den drei Parametern abhängen. In einigen Fällen sind sie auch nicht vorhanden. Peter Flom repariert den df, aber dies kann auch als unbekannter Parameter angesehen werden.
1
@PeterFlom Im Fall der Cauchy-Verteilung ist explizit angegeben, dass m und s der Ort und die Skala sind. Ich stimme der Notation m und s zu, dass sie den Mittelwert bzw. die Standardabweichung darstellen. Aber dies kann nur eine Vereinfachung \muund \sigmaals gut. +1 übrigens vor langer Zeit.
1
@PeterFlom Bedeutet dieses Zitat aus der R-Hilfedatei, dass df für die Verteilung der Schüler immer 9 ist? Denken Sie nicht, dass df auch geschätzt werden sollte? Tatsächlich ist das Fehlen von dfdie Ursache für den Fehler, und die richtige Antwort sollte ein Rezept liefern, um ihn zu finden.
Sergey Bushmanov
1
@ PeterFlom Übrigens, wenn Sie die Hilfedatei ein paar Zeilen über Ihrem Zitat lesen, werden Sie feststellen, warum sie df=9in ihrem Beispiel gut und irrelevant ist.
Sergey Bushmanov