R: nls2 错过了解决方案
R: nls2 misses the solution
我正在尝试拟合双指数函数:
t = seq(0, 30, by = 0.1)
A = 20 ; B = 10 ; alpha = 0.25 ; beta = 0.01
y = A*exp(-alpha*t) + B*exp(-beta*(t))
df = as.data.frame(cbind(t,y))
ggplot(df, aes(t, y)) + geom_line() + scale_y_continuous(limits=c(0, 50))
这个问题无法通过像log这样的简单转换来解决,所以我想使用nls2
包:
library(nls2)
fo <- y ~ Ahat*exp(-alphahat*t) + Bhat*exp(-betahat*t)
fit <- nls2(fo,
start = list(Ahat=5, Bhat=5, alphahat=0.5,betahat=0.5),
algorithm = "brute-force",
trace = TRUE,
lower = c(Ahat=0, Bhat=0, alphahat=0, betahat=0),
upper = c(Ahat=50, Bhat=50, alphahat=10,betahat=10))
fit
结果如下:
Nonlinear regression model
model: y ~ Ahat * exp(-alphahat * t) + Bhat * exp(-betahat * t)
data: NULL
Ahat Bhat alphahat betahat
5.0 5.0 0.5 0.5
residual sum-of-squares: 37910
Number of iterations to convergence: 4
Achieved convergence tolerance: NA
我假设我的代码有问题,因为:
- 数据:空?
- 为什么只有 4 次迭代?
- 很难想象 nls2 没有找到比起点更好的解决方案。
- 结果与解决方案相去甚远
您的数据是 null
因为您没有在 nls2
语句中添加任何数据。
nls2
需要这样设置:
nls2(formula, data = parent.frame(), start, control = nls.control(),
algorithm = c("default", "plinear", "port", "brute-force",
"grid-search", "random-search", "plinear-brute", "plinear-random"),
trace = FALSE, weights, ..., all = FALSE)
查看 the official documentation 以获得完整示例。
根据文档,start
参数应该是 data.frame
定义要搜索的网格的两行,或者 data.frame
包含与参数组合相对应的更多行测试您是否正在使用蛮力。此外,nls
会出现问题,因为它是完美的曲线,没有噪音。蛮力方法很慢,所以这里是一个示例,其中搜索 space 减少了 nls2
。在向数据添加一点点噪声之后,暴力 nls2
的结果然后被用作 nls
默认算法的起始值(或者您可以使用 nls2
)。
## Data
t = seq(0, 30, by = 0.1)
A = 20 ; B = 10 ; alpha = 0.25 ; beta = 0.01
y = A*exp(-alpha*t) + B*exp(-beta*(t))
df = as.data.frame(cbind(t,y))
library(nls2)
fo <- y ~ Ahat*exp(-alphahat*t) + Bhat*exp(-betahat*t)
## Define the grid to search in,
## Note: decreased the grid size
grd <- data.frame(Ahat=c(10,30),
Bhat=c(10, 30),
alphahat=c(0,2),
betahat=c(0,1))
## Do the brute-force
fit <- nls2(fo,
data=df,
start = grd,
algorithm = "brute-force",
control=list(maxiter=100))
coef(fit)
# Ahat Bhat alphahat betahat
# 10.0000000 23.3333333 0.0000000 0.3333333
## Now, run through nls:
## Fails, because there is no noise
final <- nls(fo, data=df, start=as.list(coef(fit)))
## Add a little bit of noise
df$y <- df$y+rnorm(nrow(df),0,0.001)
coef((final <- nls(fo, data=df, start=as.list(coef(fit)))))
# Ahat Bhat alphahat betahat
# 10.00034000 19.99956016 0.01000137 0.25000966
## Plot
plot(df, col="steelblue", pch=16)
points(df$t, predict(final), col="salmon", type="l")
我正在尝试拟合双指数函数:
t = seq(0, 30, by = 0.1)
A = 20 ; B = 10 ; alpha = 0.25 ; beta = 0.01
y = A*exp(-alpha*t) + B*exp(-beta*(t))
df = as.data.frame(cbind(t,y))
ggplot(df, aes(t, y)) + geom_line() + scale_y_continuous(limits=c(0, 50))
这个问题无法通过像log这样的简单转换来解决,所以我想使用nls2
包:
library(nls2)
fo <- y ~ Ahat*exp(-alphahat*t) + Bhat*exp(-betahat*t)
fit <- nls2(fo,
start = list(Ahat=5, Bhat=5, alphahat=0.5,betahat=0.5),
algorithm = "brute-force",
trace = TRUE,
lower = c(Ahat=0, Bhat=0, alphahat=0, betahat=0),
upper = c(Ahat=50, Bhat=50, alphahat=10,betahat=10))
fit
结果如下:
Nonlinear regression model
model: y ~ Ahat * exp(-alphahat * t) + Bhat * exp(-betahat * t)
data: NULL
Ahat Bhat alphahat betahat
5.0 5.0 0.5 0.5
residual sum-of-squares: 37910
Number of iterations to convergence: 4
Achieved convergence tolerance: NA
我假设我的代码有问题,因为:
- 数据:空?
- 为什么只有 4 次迭代?
- 很难想象 nls2 没有找到比起点更好的解决方案。
- 结果与解决方案相去甚远
您的数据是 null
因为您没有在 nls2
语句中添加任何数据。
nls2
需要这样设置:
nls2(formula, data = parent.frame(), start, control = nls.control(),
algorithm = c("default", "plinear", "port", "brute-force",
"grid-search", "random-search", "plinear-brute", "plinear-random"),
trace = FALSE, weights, ..., all = FALSE)
查看 the official documentation 以获得完整示例。
根据文档,start
参数应该是 data.frame
定义要搜索的网格的两行,或者 data.frame
包含与参数组合相对应的更多行测试您是否正在使用蛮力。此外,nls
会出现问题,因为它是完美的曲线,没有噪音。蛮力方法很慢,所以这里是一个示例,其中搜索 space 减少了 nls2
。在向数据添加一点点噪声之后,暴力 nls2
的结果然后被用作 nls
默认算法的起始值(或者您可以使用 nls2
)。
## Data
t = seq(0, 30, by = 0.1)
A = 20 ; B = 10 ; alpha = 0.25 ; beta = 0.01
y = A*exp(-alpha*t) + B*exp(-beta*(t))
df = as.data.frame(cbind(t,y))
library(nls2)
fo <- y ~ Ahat*exp(-alphahat*t) + Bhat*exp(-betahat*t)
## Define the grid to search in,
## Note: decreased the grid size
grd <- data.frame(Ahat=c(10,30),
Bhat=c(10, 30),
alphahat=c(0,2),
betahat=c(0,1))
## Do the brute-force
fit <- nls2(fo,
data=df,
start = grd,
algorithm = "brute-force",
control=list(maxiter=100))
coef(fit)
# Ahat Bhat alphahat betahat
# 10.0000000 23.3333333 0.0000000 0.3333333
## Now, run through nls:
## Fails, because there is no noise
final <- nls(fo, data=df, start=as.list(coef(fit)))
## Add a little bit of noise
df$y <- df$y+rnorm(nrow(df),0,0.001)
coef((final <- nls(fo, data=df, start=as.list(coef(fit)))))
# Ahat Bhat alphahat betahat
# 10.00034000 19.99956016 0.01000137 0.25000966
## Plot
plot(df, col="steelblue", pch=16)
points(df$t, predict(final), col="salmon", type="l")