[R] Fitdistr and likelihood

[Martin Maechler]

       {BCC'ed to VR's maintainer}

>>>>> "Carsten" == Carsten Steinhoff <carsten.steinhoff at stud.uni-goettingen.de>
>>>>>     on Tue, 5 Apr 2005 17:31:04 +0200 writes:

    Carsten> Hi all, I'm using the function "fitdistr" (library
    Carsten> MASS) to fit a distribution to given data.  What I
    Carsten> have to do further, is getting the
    Carsten> log-Likelihood-Value from this estimation.

    Carsten> Is there any simple possibility to realize it?

yes.  Actually, internally in fitdistr(), everything's already there.
So you need to modify fitdistr() only slightly to also return
the log likelihood.  Furthermore, subsequent implementation of a
logLik.fitdistr() method is trivial.

Here are is the unified diff against VR_7.2-14
(VR/MASS/R/fitdistr.R) :

--- /u/maechler/R/MM/STATISTICS/fitdistr.R	2004-01-22 02:16:04.000000000 +0100
+++ /u/maechler/R/MM/STATISTICS/fitdistr-improved.R	2005-04-06 10:20:25.000000000 +0200
@@ -1,3 +1,5 @@
+#### This is from VR_7.2-14  VR/MASS/R/fitdistr.R  [improved by MM]
+
 fitdistr <- function(x, densfun, start, ...)
 {
     myfn <- function(parm, ...) -sum(log(dens(parm, ...)))
@@ -11,6 +13,7 @@
         stop("'x' must be a non-empty numeric vector")
     if(missing(densfun) || !(is.function(densfun) || is.character(densfun)))
         stop("'densfun' must be supplied as a function or name")
+    n <- length(x)
     if(is.character(densfun)) {
         distname <- tolower(densfun)
         densfun <-
@@ -34,12 +37,13 @@
         if(distname == "normal") {
             if(!is.null(start))
                 stop("supplying pars for the Normal is not supported")
-            n <- length(x)
             sd0 <- sqrt((n-1)/n)*sd(x)
-            estimate <- c(mean(x), sd0)
+            mx <- mean(x)
+            estimate <- c(mx, sd0)
             sds <- c(sd0/sqrt(n), sd0/sqrt(2*n))
             names(estimate) <- names(sds) <- c("mean", "sd")
-            return(structure(list(estimate = estimate, sd = sds),
+            return(structure(list(estimate = estimate, sd = sds, n = n,
+				  loglik = sum(dnorm(x, mx, sd0, log=TRUE))),
                              class = "fitdistr"))
         }
         if(distname == "weibull" && is.null(start)) {
@@ -89,15 +93,28 @@
             parse(text = paste("densfun(x,",
                   paste("parm[", 1:l, "]", collapse = ", "),
                   ", ...)"))
+    res <-
     if("log" %in% args)
-        res <- optim(start, mylogfn, x = x, hessian = TRUE, ...)
+            optim(start, mylogfn, x = x, hessian = TRUE, ...)
     else
-        res <- optim(start, myfn, x = x, hessian = TRUE, ...)
+            optim(start, myfn, x = x, hessian = TRUE, ...)
     if(res$convergence > 0) stop("optimization failed")
     sds <- sqrt(diag(solve(res$hessian)))
-    structure(list(estimate = res$par, sd = sds), class = "fitdistr")
+    structure(list(estimate = res$par, sd = sds,
+                   loglik = - res$value, n = n), class = "fitdistr")
 }
 
+logLik.fitdistr <- function(object, REML = FALSE, ...)
+{
+    if (REML) stop("only 'REML = FALSE' is implemented")
+    val <- object$loglik
+    attr(val, "nobs") <- object$n
+    attr(val, "df") <- length(object$estimate)
+    class(val) <- "logLik"
+    val
+}
+
+
 print.fitdistr <-
     function(x, digits = getOption("digits"), ...)
 {