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