我有一个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
基本上问题是我想对分组数据进行拟合并根据拟合系数进行预测。
我nlsLM
从library(minpack.lm)
我得到一个错误
nlsModel(formula, mf, start, wts, upper) 中的错误:初始参数估计处的奇异梯度矩阵
根据@Roland,乍一看可能是模型错误或我的起始值不好。另一方面,我可以只用两个拟合参数来拟合这个模型。当我想将third
参数添加到拟合函数时,就会出现问题。
尝试 2
在那篇文章中,关注@G. Grothendieck 的建议,我尝试nlxb
从nlmrt
包中将其中一个参数固定d
为d=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”类的对象
没有!coef
或predict 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 个),阅读书籍并遵循建议。