[R] Fitting distributions
Prof Brian Ripley
ripley at stats.ox.ac.uk
Wed Sep 5 09:22:10 CEST 2001
This has come up a couple of times recently. MASS4 will have a
general-purpose function for maximum-likelihood fitting and a set of
wrappers for common distributions.
Here's a preview of the R version (rapidly converted this morning).
You can do things like
fitdist(x, dnorm, list(mean = 0, sd = 1))
fitdist(x, dnorm, list(sd = 1), mean=0.5)
the latter fixing the mean.
This is not much tested (in R and much of it had to be changed for R) but
should give a flavour of what is possible. Eventually there will be
a much more comprehensive interface.
fitdist <- function(x, dist, start, ...)
{
myfn <- function(parm, ...) -sum(log(dens(parm, ...)))
dens1 <- function(parm, x, ...) dist(x, parm, ...)
dens2 <- function(parm, x, ...) dist(x, parm[1], parm[2], ...)
dens3 <- function(parm, x, ...) dist(x, parm[1], parm[2], parm[3], ...)
dens4 <- function(parm, x, ...)
dist(x, parm[1], parm[2], parm[3], , parm[4], ...)
if(missing(x) || length(x) == 0 || mode(x) != "numeric")
stop("`x' must be a non-empty numeric vector")
if(missing(dist) || !is.function(dist))
stop("`dist' must be supplied as a function")
if(missing(start) || !is.list(start))
stop("`start' must be a named list")
nm <- names(start)
## reorder arguments to dist
f <- formals(dist)
args <- names(f)
m <- match(nm, args)
if(any(is.na(m)))
stop("`start' specifies names which are not arguments to `dist'")
formals(dist) <- c(f[c(1, m)], f[-c(1, m)])
## try to be somewhat efficient
switch(length(nm),
"1" = {dens <- dens1},
"2" = {dens <- dens2},
"3" = {dens <- dens3},
"4" = {dens <- dens3},
stop("only handle 1:4 variable parms"))
res <- optim(start, myfn, x=x, ...)
if(res$convergence > 0) stop("optimization failed")
res$par
}
--
Brian D. Ripley, ripley at stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272860 (secr)
Oxford OX1 3TG, UK Fax: +44 1865 272595
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list