-1

我正在使用用 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
}
4

1 回答 1

1

没有公式。nls正在重复传递对正在拉入“currModel”和参数（可能是 theta 和 d）的函数的调用，这些参数将根据从rescomp函数返回的标量最小化。您抱怨键入“rescomp”的结果“写在这里太长”只是表明a）您不明白您应该编辑您的问题而不是在评论中写出该输出，并且b）你对正在发生的事情的期望太窄了。

为了说明将您的第一个 nls 问题编写为函数形式：

myfunc <- function(Const,A,B,y=y,x=x) { abs(y - ( Const + A * exp(B * x)))}

因此，您将最小化 y 与预测的绝对偏差：

nlmod <- nls( ~myfunc(Const,A,B,y=y,x=x) ,start=list(Const=5,A=5,B=7))
nlmod  # same results
于 2013-10-09T18:59:33.200 回答