Estymator NW w R
W statystyce matematycznej bardzo często stosowany jest estymator największej wiarygodności, którego idea polega na poszukiwaniu maksimum tzw. funkcji wiarygodności.
W artykule omówiona zostanie biblioteka maxLik, w której zaimplementowane są różne algorytmy numeryczne poszukiwania NW-estymatora.
Opis biblioteki znajdujesię na stronie http://cran.r-project.org/web/packages/maxLik/index.html
W pierwszej kolejności należy bibliotekę zainstalować. Można to zrobić np. za pomocą polecenia
install.packages( "maxLik", repos="http://R-Forge.R-project.org" )
gdzie repos to adres miejsca z którego ma być pobrana biblioteka.
Przed użyciem po raz pierwszy bibliotekę trzeba załadować. Robimy to poleceniem
library("maxLik")
Do poszukiwania NW-estymatora aktualnie możemy użyć następujących algorytmów
- Newton-Raphson (maxNR)
- Berndt-Hall-Hall-Hausman (maxBHHH)
- Broyden-Fletcher-Goldfarb-Shanno (maxBFGS)
- simulated annealing (maxSANN)
- Nelder-Mead (maxNM)
W nawiasach podano nazwy funkcji dostępnych w R.
Dodatkowo dostępna jest finkcja maxLik za pomocą w prosty sposób (zmieniając tylko wartość jednego argumentu można zbadać wyniki otrzymana z różnych algorytmów).Funkcji ma 5 argumentów. Pierwszy argument tej funkcji to loglik Jest on obowiązkowy i określa funkcje wiarygodności. Jego pierwszy argument (funkcji loglik) musi być wektor parametrów które mają być oszacowane. Drugim i trzecim argumentem maxLik (grad i hess - opcjonalnie) są funkcje, które zwracają gradient i hesian dla logarytmu fiunkcji wiarygodności. Jeśli te funkcje nie są podawane przez użytkownika, numeryczne gradienty i hessiany obliczane są w razie potrzeby. Czwarty argument (start) jest obowiązkowy i jest używany do zaincjowania algorytmu. Wreszcie, piąty argument (method) jest opcjonalny i pozwala na wybór jednego z algorytmów: "NR" (domyślna wartośc), "BHHH", "BFGS", "NM" lub "SANN".
Obiekt zwracany przez funkcję maxLik składa się z:
- activePar -wartość logiczna wskazujaca czy parametrbył staly
- estimate - estymator NW
- hessian - wartość hessianu dla wyznaczonego NW-estymatora
- last.step - lista przedstawiająca wartości dla ostatniego kroku (w przypadku funkcji nieograniczonych)
- message - opis otrzmanego wyniku
- code -liczba calkowita symbolizująca warunek stopu, który zastał użyty
- gradient - wartość gradientu dla NW-estymatora
- iterations - ilość rzeczywiście wykonanych iteracji
- maximum - wartość logarytmu wiarygodności dla NW-estymatora
- type - opis użytej metody optymalizacji
Przykład 1
## estymacja parametru skali (theta) rozkładu wykladniczego
t <- rexp(100, 2) #losowniepróby
loglik <- function(theta) log(theta) - theta*t #logarytm wiarygodności
gradlik <- function(theta) 1/theta - t # gradient logarytmu wiarygodności
hesslik <- function(theta) -100/theta^2 # hesian logarytmu wiarygodności
## estymacja z punktem startowym 1
a <- maxLik(loglik, start=1, print.level=2, method="NR")
a
a$estimate
## estymacja z gradientem i hessianem podanym w sposób analityczny
a <- maxLik(loglik, gradlik, hesslik, start=1)
a
Przykład 2
## Przykład z argumentem wektorowym: estymacja średniej i wariancji rozkładu normalnego ##za pomocą NW-estymatora
loglik <- function(param) {
mu <- param[1]
sigma <- param[2]
ll <- -0.5*N*log(2*pi) - N*log(sigma) - sum(0.5*(x - mu)^2/sigma^2)
ll
}
x <- rnorm(1000, 1, 2) # use mean=1, stdd=2
N <- length(x)
res <- maxLik(loglik, start=c(1,1))
res
Załącznik | Wielkość |
---|---|
maxLik.r | 1015 bajtów |