全部版块 我的主页
论坛 数据科学与人工智能 数据分析与数据科学 R语言论坛
3276 4
2011-03-07
我调用了程序包里的一个函数hyperbFit()拟合dataVector数据的分布,结果出来了,如下:

Data:      dataVector
Parameter estimates:
       mu      delta      alpha       beta  
0.001160   0.001918  86.117240  -3.427075  
Likelihood:         8201.185
Method:             Nelder-Mead
Convergence code:   0
Iterations:         375

其中,mu      delta      alpha       beta  就是拟合的分布参数,我想提取出这几个参数用于继续编程,可是内存里没有这几个对象,怎么才能提取这4个参数呢?

下面是这个函数的代码,高手帮我分析一下,难道别人包里的数据无法提取,就像C语言函数内部的变量一样?

> hyperbFit

function (x, freq = NULL, breaks = NULL, paramStart = NULL, startMethod = "Nelder-Mead",
    startValues = "BN", method = "Nelder-Mead", hessian = FALSE,
    plots = FALSE, printOut = FALSE, controlBFGS = list(maxit = 200),
    controlNM = list(maxit = 1000), maxitNLM = 1500, ...)
{
    xName <- paste(deparse(substitute(x), 500), collapse = "\n")
    if (!is.null(freq)) {
        if (length(freq) != length(x))
            stop("vectors x and freq are not of the same length")
        x <- rep(x, freq)
    }
    x <- as.numeric(na.omit(x))
    startInfo <- hyperbFitStart(x, breaks = breaks, startValues = startValues,
        paramStart = paramStart, startMethodSL = startMethod,
        startMethodMoM = startMethod, ...)
    paramStart <- startInfo$paramStart
    svName <- startInfo$svName
    breaks <- startInfo$breaks
    empDens <- startInfo$empDens
    midpoints <- startInfo$midpoints
    llfunc <- function(param) {
        mu <- param[1]
        delta <- param[2]
        alpha <- param[3]
        beta <- param[4]
        hyperbPi <- beta/sqrt(alpha^2 - beta^2)
        zeta <- delta * sqrt(alpha^2 - beta^2)
        KNu <- besselK(zeta, nu = 1)
        hyperbDens <- (2 * delta * sqrt(1 + hyperbPi^2) * KNu)^(-1) *
            exp(-zeta * (sqrt(1 + hyperbPi^2) * sqrt(1 + ((x -
                mu)/delta)^2) - hyperbPi * (x - mu)/delta))
        as.numeric(hyperbDens)
        -sum(log(hyperbDens))
    }
    output <- numeric(7)
    ind <- 1:4
    if (method == "BFGS") {
        opOut <- optim(paramStart, llfunc, NULL, method = "BFGS",
            control = controlBFGS, ...)
    }
    if (method == "Nelder-Mead") {
        opOut <- optim(paramStart, llfunc, NULL, method = "Nelder-Mead",
            control = controlNM, ...)
    }
    if (method == "nlm") {
        ind <- c(2, 1, 5, 4)
        opOut <- nlm(llfunc, paramStart, iterlim = maxitNLM,
            ...)
    }
    param <- as.numeric(opOut[[ind[1]]])[1:4]
    names(param) <- c("mu", "delta", "alpha", "beta")
    maxLik <- -(as.numeric(opOut[[ind[2]]]))
    conv <- as.numeric(opOut[[ind[4]]])
    iter <- as.numeric(opOut[[ind[3]]])[1]
    if (hessian) {
        n <- length(param)
        lloutput <- llfunc(param)
        eps <- .Machine$double.eps
        h = eps^(1/3) * apply(as.data.frame(param), 1, function(z) max(abs(z),
            0.01))
        ee = diag(h)
        gp = vector(mode = "numeric", length = n)
        gm = vector(mode = "numeric", length = n)
        for (i in 1:n) gp[i] <- llfunc(param + ee[, i])
        for (i in 1:n) gm[i] <- llfunc(param - ee[, i])
        H = h %*% t(h)
        Hm = H
        Hp = H
        for (i in 1:n) {
            for (j in i:n) {
                Hp[i, j] <- llfunc(param + ee[, i] + ee[, j])
                Hp[j, i] <- Hp[i, j]
                Hm[i, j] <- llfunc(param - ee[, i] - ee[, j])
                Hm[j, i] <- Hm[i, j]
            }
        }
        for (i in 1:n) {
            for (j in i:n) {
                H[i, j] = ((Hp[i, j] - gp[i] - gp[j] + lloutput +
                  lloutput - gm[i] - gm[j] + Hm[i, j])/H[i, j])/2
                H[j, i] = H[i, j]
            }
        }
        colnames(H) <- names(param)
        rownames(H) <- names(param)
        opOut$hessian <- H
    }
    fitResults <- list(param = param, maxLik = maxLik, hessian = if (hessian) opOut$hessian else NULL,
        method = method, conv = conv, iter = iter, obs = x, obsName = xName,
        paramStart = paramStart, svName = svName, startValues = startValues,
        breaks = breaks, midpoints = midpoints, empDens = empDens)
    class(fitResults) <- c("hyperbFit", "distFit")
    if (printOut)
        print(fitResults, ...)
    if (plots)
        plot.hyperbFit(fitResults, ...)
    fitResults
}
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

全部回复
2011-3-8 01:09:27
假设你调用函数的返回值是myfit
用myfit$param就得到了返回的参数值
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2011-3-9 23:09:31
多谢 2# qoiqpwqr
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2014-9-8 22:34:39
qoiqpwqr 发表于 2011-3-8 01:09
假设你调用函数的返回值是myfit
用myfit$param就得到了返回的参数值
刚刚用您的回复解决了困扰我的问题,太感谢了~~~
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2014-9-8 22:34:40
qoiqpwqr 发表于 2011-3-8 01:09
假设你调用函数的返回值是myfit
用myfit$param就得到了返回的参数值
刚刚用您的回复解决了困扰我的问题,太感谢了~~~
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

相关推荐
栏目导航
热门文章
推荐文章

说点什么

分享

扫码加好友,拉您进群
各岗位、行业、专业交流群