发帖
雷达卡
我调用了程序包里的一个函数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
}
京公网安备 11010802022788号
