Pas t-distribusi di R: parameter penskalaan

Bagaimana saya menyesuaikan parameter distribusi-t, yaitu parameter yang sesuai dengan 'rata-rata' dan 'standar deviasi' dari distribusi normal. Saya menganggap mereka disebut 'berarti' dan 'scaling / derajat kebebasan' untuk distribusi-t?

Kode berikut sering menghasilkan kesalahan 'optimasi gagal'.

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

Apakah saya harus mengukur x terlebih dahulu atau mengubahnya menjadi probabilitas? Bagaimana cara terbaik untuk melakukannya?

fitdistrmenggunakan kemungkinan-maksimum dan teknik optimasi untuk menemukan parameter dari distribusi yang diberikan. Terkadang, terutama untuk distribusi-t, seperti yang diperhatikan oleh @ user12719, optimasi dalam bentuk:

fitdistr(x, "t")

gagal dengan kesalahan.

Dalam hal ini Anda harus memberikan pengoptimal tangan dengan memberikan titik awal dan batas bawah untuk mulai mencari parameter yang optimal:

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

Catatan, df=3 adalah tebakan terbaik Anda untuk apa yang "optimal" dfbisa. Setelah memberikan info tambahan ini kesalahan Anda akan hilang.

Beberapa kutipan untuk membantu Anda lebih memahami mekanisme dalam fitdistr :

Untuk distribusi Normal, log-Normal, geometris, eksponensial, dan Poisson, MLEs bentuk-tertutup (dan kesalahan standar yang tepat) digunakan, dan start tidak boleh disediakan.

...

Untuk distribusi bernama berikut, nilai awal yang masuk akal akan dihitung jika startdihilangkan atau hanya ditentukan sebagian: "cauchy", "gamma", "logistic", "binomial negatif" (ditentukan oleh ukuran dan ukuran), "t" dan "weibull ". Perhatikan bahwa nilai awal ini mungkin tidak cukup baik jika kecocokannya buruk: khususnya mereka tidak tahan terhadap pencilan kecuali distribusi yang dipasang berekor panjang.

$\nu$$t$ , 1 ).

$\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$ besar. Tetapi apakah estimator kemungkinan maksimum itu ada gunanya?

Mari kita coba beberapa simulasi:

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

Menunjukkan estimasi sangat tidak stabil (melihat histogram, sebagian besar nilai estimasi berada di batas atas yang diberikan untuk mengoptimalkan 200).

Mengulangi dengan ukuran sampel yang lebih besar:

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

yang jauh lebih baik, tetapi rata-rata masih jauh di atas nilai sebenarnya 2,5.

Kemudian ingat bahwa ini adalah versi sederhana dari masalah sebenarnya di mana parameter lokasi dan skala juga harus diperkirakan.

$t$$\nu$

Dalam bantuan untuk fitdistr adalah contoh ini:

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

menunjukkan bahwa Anda hanya perlu nilai untuk df. Tapi itu mengasumsikan standardisasi.

Untuk kontrol yang lebih besar, mereka juga tampil

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

di mana parameternya adalah m = rata-rata, s = standar deviasi, df = derajat kebebasan