加载中…
个人资料
  • 博客等级:
  • 博客积分:
  • 博客访问:
  • 关注人气:
  • 获赠金笔:0支
  • 赠出金笔:0支
  • 荣誉徽章:
正文 字体大小:

R软件中的AIC函数

(2012-03-15 17:35:10)
标签:

杂谈

分类: R软件学习

#  File src/library/stats/R/AIC.R
#  Part of the R package, http://www.R-project.org
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation; either version 2 of the License, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
#  A copy of the GNU General Public License is available at
#  http://www.r-project.org/Licenses/

#### Return the value of Akaike's Information Criterion

AIC <- function(object, ..., k = 2) UseMethod("AIC")

## For back-compatibility
AIC.logLik <- function(object, ..., k = 2)
    -2 * as.numeric(object) + k * attr(object, "df")

AIC.default <- function(object, ..., k = 2)
{
    ## AIC for various fitted objects --- any for which there's a logLik() method:
    ll <- if("stats4" %in% loadedNamespaces()) stats4:::logLik else logLik
    if(length(list(...))) {# several objects: produce data.frame
        lls <- lapply(list(object, ...), ll)
        vals <- sapply(lls, function(el) {
            no <- attr(el, "nobs")
            c(as.numeric(el), attr(el, "df"),
              if(is.null(no)) NA_integer_ else no)
        })
        val <- data.frame(df = vals[2L,], ll = vals[1L,])
        nos <- na.omit(vals[3L,])
        if (length(nos) && any(nos != nos[1L]))
            warning("models are not all fitted to the same number of observations")
        val <- data.frame(df = val$df, AIC = -2*val$ll + k*val$df)
        Call <- match.call()
        Call$k <- NULL
        row.names(val) <- as.character(Call[-1L])
        val
    } else {
        lls <- ll(object)
         -2 * as.numeric(lls) + k * attr(lls, "df")
    }
}

BIC <- function(object, ...) UseMethod("BIC")

## For back-compatibility
BIC.logLik <- function(object, ...)
    -2 * as.numeric(object) + attr(object, "df") * log(nobs(object))

BIC.default <- function(object, ...)
{
    ll <- if("stats4" %in% loadedNamespaces()) stats4:::logLik else logLik
    Nobs <- if("stats4" %in% loadedNamespaces()) stats4:::nobs else nobs
    if(length(list(...))) {
        lls <- lapply(list(object, ...), ll)
        vals <- sapply(lls, function(el) {
            no <- attr(el, "nobs")
            c(as.numeric(el), attr(el, "df"),
              if(is.null(no)) NA_integer_ else no)
        })
        val <- data.frame(df = vals[2L,], ll = vals[1L,], nobs = vals[3L,])
        nos <- na.omit(val$nobs)
        if (length(nos) && any(nos != nos[1L]))
            warning("models are not all fitted to the same number of observations")
        ## if any val$nobs = NA, try to get value via nobs().
        unknown <- is.na(val$nobs)
        if(any(unknown))
            val$nobs[unknown] <-
                sapply(list(object, ...)[unknown],
                       function(x, f) tryCatch(f(x), error = function(e) NA_real_),
                       f = Nobs)
        val <- data.frame(df = val$df, BIC = -2*val$ll + log(val$nobs)*val$df)
        row.names(val) <- as.character(match.call()[-1L])
        val
    } else {
        lls <- ll(object)
        nos <- attr(lls, "nobs")
        if (is.null(nos))
            nos <- tryCatch(Nobs(object), error = function(e) NA_real_)
        -2 * as.numeric(lls) + log(nos) * attr(lls, "df")
    }
}

0

阅读 收藏 喜欢 打印举报/Report
  

新浪BLOG意见反馈留言板 欢迎批评指正

新浪简介 | About Sina | 广告服务 | 联系我们 | 招聘信息 | 网站律师 | SINA English | 产品答疑

新浪公司 版权所有