【问题标题】:nls function in R within functions and listsR中的nls函数在函数和列表中
【发布时间】: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


【解决方案1】:

没有公式。 nls 反复调用一个函数,该函数正在拉入 'currModel' 和参数(可能是 theta 和 d),这些参数将根据从 rescomp 函数返回的标量最小化。您抱怨键入“rescomp”的结果“写在这里太长”只是表明a)您不明白您应该编辑您的问题而不是在 cmets 中写入该输出,并且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

【讨论】:

  • 感谢您的回答,真的很有帮助。我还在上面添加了rescomp函数
  • 如果您想查看驱动模型构建的公式,那么我建议您在rescomp 中找不到它,而是通过剖析currModel 对象。我的第一个建议是尝试currModel$call ...如果不成功,请报告names(currModel)class(currModel) 的结果。
猜你喜欢
  • 2013-02-08
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2022-01-09
  • 1970-01-01
  • 2013-12-05
  • 2021-11-20
  • 1970-01-01
相关资源
最近更新 更多