{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"), ...)
{