无法使用非线性拟合方法（nlsLM、nlxb 和 wrapnls）进行拟合
Failing to do fitting with non linear fitting methods (nlsLM, nlxb and wrapnls)
我有一个 nls fitting 任务想用 R 完成。
我第一次尝试这样做 here 正如@Roland 指出的
"The point is that complex models are difficult to fit. The more so, the less the data supports the model until it become impossible. You might be able to fit this, if you had extremely good starting values."
我同意@Roland 但如果 excel 可以做到这一点为什么 R 不能做到？
基本上这个拟合可以用Excel的GRG Nonlinear solver来完成，但是这个过程非常耗时，而且有时拟合不好。 （因为现实中有很多数据）
这是我的样本 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
我已经 post 在这里提问 trouble-when-adding-3rd-fitting-parameter-in-nls
基本上问题是我想对分组数据进行拟合并根据拟合系数进行预测。
我使用了 library(minpack.lm) 中的 nlsLM 我得到了一个错误
Error in nlsModel(formula, mf, start, wts, upper) :
singular gradient matrix at initial parameter estimates
根据@Roland 的说法，乍一看可能是模型错误或我的起始值不好。另一方面，我可以只用两个拟合参数来拟合这个模型。当我想将 third 参数添加到拟合函数时，问题就出现了。
尝试 2
其中 post trouble-when-adding-3rd-fitting-parameter-in-nls 通过关注@G。 Grothendieck 建议，我尝试了 nlmrt 包中的 nlxb 并将参数之一 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
拟合系数合理 wolaa！但是这次我意识到（@G.Grothendieck 后来也指出了）nlxb 之后不可能预测新的值（为什么=？我不知道！）
predvals <- ldply(fit, .fun=predictvals, xvar="time", yvar="value",xrange=range(range)) # predict values
::您可以从 here
中找到 predictvals 函数
Error in UseMethod("predict") :
no applicable method for 'predict' applied to an object of class "nlmrt"
没有！ coef 或 predict methods 用于 "nlmrt" class 个对象。
尝试 3
关注@G后。格洛腾迪克的另一个建议
接下来我尝试从 nlmrt.
wrapnls
因为在这篇post中他说，
can-we-make-prediction-with-nlxb-from-nlmrt-package
”因为 nlmrt 包确实提供了 wrapnls 它将 运行 nlmrt 然后是 nls 这样一个 "nls" 对象结果然后该对象可以与所有 "nls" class 方法一起使用。
来自同一个 nlmrt 包仍然有如下问题
我在第一个 post 后放弃使用 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))
Error in nlsModel(formula, mf, start, wts, upper) :
singular gradient matrix at initial parameter estimates
我回到了开始的地方，出现了这个错误。
我想我已经尽我所能，寻找相关示例（尽管只有 3 个），阅读书籍并遵循建议。
像这样在 nlxb 之后使用 nls2 包中的 nls2（假设 data_rep，
formula、d_step、f 和 d 来自问题）。为了使示例最小化，我们删除了 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" class 对象，其系数与 fit.nlxb 相同，我们可以使用 fitted() 和 predict() 以及所有其他 "nls" 方法。
我有一个 nls fitting 任务想用 R 完成。
我第一次尝试这样做 here 正如@Roland 指出的
"The point is that complex models are difficult to fit. The more so, the less the data supports the model until it become impossible. You might be able to fit this, if you had extremely good starting values."
我同意@Roland 但如果 excel 可以做到这一点为什么 R 不能做到？
基本上这个拟合可以用Excel的GRG Nonlinear solver来完成，但是这个过程非常耗时，而且有时拟合不好。 （因为现实中有很多数据）
这是我的样本 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
我已经 post 在这里提问 trouble-when-adding-3rd-fitting-parameter-in-nls
基本上问题是我想对分组数据进行拟合并根据拟合系数进行预测。
我使用了 library(minpack.lm) 中的 nlsLM 我得到了一个错误
Error in nlsModel(formula, mf, start, wts, upper) : singular gradient matrix at initial parameter estimates
根据@Roland 的说法，乍一看可能是模型错误或我的起始值不好。另一方面，我可以只用两个拟合参数来拟合这个模型。当我想将 third 参数添加到拟合函数时，问题就出现了。
尝试 2
其中 post trouble-when-adding-3rd-fitting-parameter-in-nls 通过关注@G。 Grothendieck 建议，我尝试了 nlmrt 包中的 nlxb 并将参数之一 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
拟合系数合理 wolaa！但是这次我意识到（@G.Grothendieck 后来也指出了）nlxb 之后不可能预测新的值（为什么=？我不知道！）
predvals <- ldply(fit, .fun=predictvals, xvar="time", yvar="value",xrange=range(range)) # predict values
::您可以从 here
中找到predictvals 函数
Error in UseMethod("predict") : no applicable method for 'predict' applied to an object of class "nlmrt"
没有！ coef 或 predict methods 用于 "nlmrt" class 个对象。
尝试 3
关注@G后。格洛腾迪克的另一个建议
接下来我尝试从 nlmrt.
wrapnls
因为在这篇post中他说， can-we-make-prediction-with-nlxb-from-nlmrt-package
”因为 nlmrt 包确实提供了 wrapnls 它将 运行 nlmrt 然后是 nls 这样一个 "nls" 对象结果然后该对象可以与所有 "nls" class 方法一起使用。
来自同一个 nlmrt 包仍然有如下问题
我在第一个 post 后放弃使用 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))
Error in nlsModel(formula, mf, start, wts, upper) : singular gradient matrix at initial parameter estimates
我回到了开始的地方，出现了这个错误。 我想我已经尽我所能，寻找相关示例（尽管只有 3 个），阅读书籍并遵循建议。
像这样在 nlxb 之后使用 nls2 包中的 nls2（假设 data_rep，
formula、d_step、f 和 d 来自问题）。为了使示例最小化，我们删除了 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" class 对象，其系数与 fit.nlxb 相同，我们可以使用 fitted() 和 predict() 以及所有其他 "nls" 方法。