【发布时间】:2013-10-09 17:13:06
【问题描述】:
我正在使用一个用 R 编写的包,以将实验数据拟合到特定模型 我想编写自己的模型并使用相同的包(作者声称这是可能的)。
所以我正在挖掘源文件以找到它们定义模型函数的位置。我被困住了
顺便说一句,他们使用nls。我了解nls 在下面的案例中是如何工作的
x <- -(1:100)/10
y <- 100 + 10 * exp(x / 2) + rnorm(x)/2
nlmod <- nls(y ~ Const + A * exp(B * x),start=list(Const=5,A=5,B=7))
创建一个数据集并在它旁边编写模型。
但是他们使用nls,如下所示:
> nls(~rescomp(theta = t, d = d, currModel = currModel),
> data = list(d = vector(), currModel = currModel))
当我跑步时
> ~rescomp(theta = t, d = d, currModel = currModel)
查看上面示例中的公式 (Const + A * exp(B * x))。我明白了
~rescomp(theta = t, d = d, currModel = currModel)
environment: 0x000000000a7f7530
我想了解如何查看结果公式以及nls
在数据设置为 list\environment 时工作。有什么建议吗?
这里是rescomp的内容
function (theta = vector(), d = vector(), currModel = currModel,
currTheta = vector())
{
if (length(currTheta) == 0)
currTheta <- getThetaCl(theta, currModel)
groups <- currModel@groups
m <- currModel@modellist
resid <- clpindepX <- list()
nexp <- length(m)
for (i in 1:nexp) {
clpindepX[[i]] <- if (!m[[i]]@clpdep || m[[i]]@getX)
getClpindepX(model = m[[i]], theta = currTheta[[i]],
multimodel = currModel, returnX = FALSE, rawtheta = theta,
dind = 0)
else matrix()
}
for (i in 1:length(groups)) {
resid[[i]] <- residPart(model = m[[1]], group = groups[[i]],
multimodel = currModel, thetalist = currTheta, clpindepX = clpindepX,
finished = currModel@finished, returnX = FALSE, rawtheta = theta)
if (currModel@finished) {
currModel <- fillResult(group = groups[[i]], multimodel = currModel,
thetalist = currTheta, clpindepX = clpindepX,
rlist = resid[[i]], rawtheta = theta)
}
}
if (currModel@finished) {
currModel@fit@nlsres$onls$nclp <- currModel@nclp
if (currModel@optlist[[1]]@sumnls) {
if (class(currModel@fit@nlsres$onls) == "nls")
class(currModel@fit@nlsres$onls) <- "timp.nls"
else if (class(currModel@fit@nlsres$onls) == "nls.lm")
class(currModel@fit@nlsres$onls) <- "timp.nls.lm"
else class(currModel@fit@nlsres$onls) <- "timp.optim"
currModel@fit@nlsres$sumonls <- summary(currModel@fit@nlsres$onls,
currModel = currModel, currTheta = currTheta)
}
if (currModel@stderrclp) {
for (i in 1:length(groups)) {
currModel <- getStdErrClp(group = groups[[i]],
multimodel = currModel, thetalist = currTheta,
clpindepX = clpindepX, rlist = resid[[i]],
rawtheta = theta)
}
if (currModel@stderrclp) {
for (i in 1:length(groups)) {
currModel <- getStdErrClp(group = groups[[i]],
multimodel = currModel, thetalist = currTheta,
clpindepX = clpindepX, rlist = resid[[i]],
rawtheta = theta)
}
}
}
if (currModel@finished && currModel@trilinear) {
trires <- triResolve(currModel, currTheta)
currModel <- trires$currModel
currTheta <- trires$currTheta
}
if (currModel@finished && m[[1]]@mod_type == "kin") {
if (m[[1]]@fullk) {
for (i in 1:nexp) {
nocolsums <- length(m[[1]]@lightregimespec) >
0
eig <- fullKF(currTheta[[i]]@kinpar, currTheta[[i]]@kinscal,
m[[1]]@kmat, currTheta[[i]]@jvec, m[[1]]@fixedkmat,
m[[1]]@kinscalspecial, m[[1]]@kinscalspecialspec,
nocolsums)
currTheta[[i]]@eigenvaluesK <- eig$values
}
}
}
if (currModel@finished) {
return(list(currModel = currModel, currTheta = currTheta))
}
if (currModel@algorithm == "optim")
retval <- sum(unlist(resid))
else retval <- unlist(resid)
retval
}
【问题讨论】:
-
在示例代码中加入
rm(list=ls())不是一个好主意。 -
尝试只输入
rescomp -
O先生,我忘记删除那行了。
-
@Ben Bolker,我已经查看了
rescomp以找到公式,但找不到。rescomp转function.rescomp的内容太长这里就不写了 -
nls只是使用不同的起始变量评估rescomp以最大限度地减少错误。
标签: r list function environment nls