1

我有一个nls fitting我想用 R 做的任务。我第一次尝试在这里做这件事,正如@Roland 指出的那样

“关键是复杂的模型很难拟合。越是这样,支持模型的数据就越少,直到它变得不可能。如果你有非常好的起始值,你也许可以拟合这个。”

我可以同意@Roland,但如果excel可以做到这一点,为什么R不能做到呢?

基本上这种拟合可以用 Excel 的 GRG 非线性求解器完成,但是这个过程非常耗时,有时拟合效果不好。(因为现实中有很多数据)。

这是我的示例data.frame。我想set用下面提供的模型来拟合每个组,

set.seed(12345)
    set =rep(rep(c("1","2","3","4"),each=21),times=1)
    time=rep(c(10,seq(100,900,100),seq(1000,10000,1000),20000),times=1)
    value <- replicate(1,c(replicate(4,sort(10^runif(21,-6,-3),decreasing=FALSE))))
    data_rep <- data.frame(time, value,set)

     > head(data_rep)
            #    time        value set
            #1     10 1.007882e-06   1
            #2    100 1.269423e-06   1
            #3    200 2.864973e-06   1
            #4    300 3.155843e-06   1
            #5    400 3.442633e-06   1
            #6    500 9.446831e-06   1
            *      *       *         *

尝试 1

我已经在这里发布了一个问题Trouble-when-adding-3rd-fitting-parameter-in-nls

基本上问题是我想对分组数据进行拟合并根据拟合系数进行预测。

nlsLMlibrary(minpack.lm)我得到一个错误

nlsModel(formula, mf, start, wts, upper) 中的错误:初始参数估计处的奇异梯度矩阵

根据@Roland,乍一看可能是模型错误或我的起始值不好。另一方面,我可以只用两个拟合参数来拟合这个模型。当我想将third参数添加到拟合函数时,就会出现问题。

尝试 2

在那篇文章,关注@G. Grothendieck 的建议,我尝试nlxbnlmrt包中将其中一个参数固定dd=32并进行如下拟合;

formula = value~Ps*(1-exp(-2*f*time*exp(-d)))*1/(sqrt(2*pi*sigma))*exp(-(d-d_ave)^2/(2*sigma))*d_step

d_step <- 1
f <- 1e9
d <- 32    
library(plyr)
library(nlmrt)
get.coefs <- function(data_rep) {
fit <-  nlxb(formula ,
              data = data_rep,
              start=c(d_ave=44,sigma=12,Ps=0.5),
              lower=c(d_ave=25,sigma=2,Ps=0.5),
              upper=c(d_ave=60,sigma=15,Ps=1),
              trace=TRUE)
}

fit <- dlply(data_rep, c("set"), .fun = get.coefs)   # Fit data grouped by "set"

#    > fit
#    $`1`
#    nlmrt class object: x 
#    residual sumsquares =  1.474e-07  on  21 observations
#        after  12    Jacobian and  13 function evaluations
#      name            coeff          SE       tstat      pval      #gradient    JSingval   
#    d_ave            42.0126            NA         NA         NA  #-7.082e-15    0.001733  
#    sigma            12.8377            NA         NA         NA   #2.408e-15   1.289e-19  
#    Ps              0.973223            NA         NA         NA    #9.33e-15    3.37e-20  
#    
#    $`2`
#    nlmrt class object: x 
#    residual sumsquares =  6.2664e-08  on  21 observations
#        after  12    Jacobian and  13 function evaluations
#      name            coeff          SE       tstat      pval      #gradient    JSingval   
#    d_ave             42.246            NA         NA         NA  #-7.269e-15    0.001428  
#    sigma            12.7429            NA         NA         NA   #2.568e-15   3.098e-19  
#    Ps              0.981517            NA         NA         NA   #9.211e-15   2.746e-20  
#    
#    $`3`
#    nlmrt class object: x 
#    residual sumsquares =  1.773e-07  on  21 observations
#        after  12    Jacobian and  13 function evaluations
#      name            coeff          SE       tstat      pval      #gradient    JSingval   
#    d_ave             41.968            NA         NA         NA  #-6.438e-15    0.001798  
#    sigma            12.8561            NA         NA         NA   #2.173e-15   2.414e-19  
#    Ps              0.972988            NA         NA         NA   #8.534e-15   5.922e-20  

#    $`4`
#    nlmrt class object: x 
#    residual sumsquares =  2.5219e-07  on  21 observations
#        after  12    Jacobian and  13 function evaluations
#      name            coeff          SE       tstat      pval      #gradient    JSingval   
#    d_ave            41.8532            NA         NA         NA  #-4.454e-15    0.001976  
#    sigma            12.9045            NA         NA         NA   #1.474e-15   3.443e-19  
#    Ps              0.974319            NA         NA         NA   #5.987e-15   3.124e-20  

#    attr(,"split_type")
#    [1] "data.frame"
#    attr(,"split_labels")
#      set
#    1   1
#    2   2
#    3   3
#    4   4

拟合系数是合理的哇!但是这一次我意识到(@G. Grothendieck 后来也指出)不可能预测之后的新值nlxb(为什么=?我不知道!)

predvals <- ldply(fit, .fun=predictvals, xvar="time", yvar="value",xrange=range(range)) # predict values 

::你可以从这里predictvals找到功能

UseMethod(“predict”)中的错误:没有适用于“predict”的方法应用于“nlmrt”类的对象

没有!coefpredict methods用于"nlmrt"类对象。

尝试 3

在关注@G 之后。格洛腾迪克接下来我尝试wrapnls了另一个建议nlmrt

因为在这篇文章中他说, can-we-make-prediction-with-nlxb-from-nlmrt-package

“因为 nlmrt 包确实提供wrapnls了哪些将运行nlmrt,然后nls生成一个"nls"对象,然后该对象可以与所有"nls"类方法一起使用。

从同一个nlmrt包中仍然遇到如下问题

plyr我在第一篇文章后放弃使用,因为加载plyr并且dplyr使我的问题更加复杂。所以我会坚持dplyr使用do函数。

library(dplyr)
library(nlmrt)
formula = value~Ps*(1-exp(-2*f*time*exp(-d)))*1/(sqrt(2*pi*sigma))*exp(-(d-d_ave)^2/(2*sigma))*d_step

d_step <- 1
f <- 1e9
d <- 32

dffit = data_rep %>% group_by(set) %>%
  do(fit = wrapnls(formula ,
                data = .,
                start=c(d_ave=44,sigma=12,Ps=0.5),
                lower=c(d_ave=25,sigma=2,Ps=0.5),
                upper=c(d_ave=60,sigma=15,Ps=1),
                trace=TRUE))

nlsModel(formula, mf, start, wts, upper) 中的错误:初始参数估计处的奇异梯度矩阵

我回到了我从这个错误开始的地方。我想我尝试了我能做的一切,寻找相关的例子(虽然只有 3 个),阅读书籍并遵循建议。

4

1 回答 1

2

nls2像这样从 nls2 包中使用nlxb(假设data_rep, formula, d_step,fd来自问题)。为了使示例最小化,我们消除了 dplyr 并仅显示 set == 2 的计算。

library(nlmrt)
library(nls2)

data_rep2 <- subset(data_rep, set == 2)

fit.nlxb <- nlxb(formula , data = data_rep2,
                start = c(d_ave = 44, sigma = 12, Ps = 0.5),
                lower = c(d_ave = 25, sigma = 2, Ps = 0.5),
                upper = c(d_ave = 60, sigma = 15, Ps = 1))

fit.nls <- nls2(formula, data = data_rep2, start = fit.nlxb$coefficients,
  algorithm = "brute-force")

identical(fit.nlxb$coefficients, coef(fit.nls))
## [1] TRUE

fit.nls是一个"nls"具有相同系数的类对象,fit.nlxb我们可以在其上使用fitted()andpredict()和所有其他"nls"方法。

于 2016-05-11T00:49:56.637 回答