强制nls拟合通过指定点的曲线
Forcing nls to fit a curve passing through a specified point
我正在尝试将 Boltzmann sigmoid 1/(1+exp((x-p1)/p2))
拟合到这个小型实验数据集:
xdata <- c(-60,-50,-40,-30,-20,-10,-0,10)
ydata <- c(0.04, 0.09, 0.38, 0.63, 0.79, 1, 0.83, 0.56)
我知道这很简单。例如,使用 nls
:
fit <-nls(ydata ~ 1/(1+exp((xdata-p1)/p2)),start=list(p1=mean(xdata),p2=-5))
我得到以下结果:
Formula: ydata ~ 1/(1 + exp((xdata - p1)/p2))
Parameters:
Estimate Std. Error t value Pr(>|t|)
p1 -33.671 4.755 -7.081 0.000398 ***
p2 -10.336 4.312 -2.397 0.053490 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1904 on 6 degrees of freedom
Number of iterations to convergence: 13
Achieved convergence tolerance: 7.079e-06
然而,我需要(由于理论原因)拟合曲线恰好通过点(-70, 0)
。虽然上面显示的拟合表达式的值在 x = -70
处接近零,但它并不完全为零,这不是我想要的。
所以,问题是:有没有办法告诉 nls
(或其他一些函数)拟合相同的表达式,但强制它通过指定的点?
更新:
正如评论中提到的那样,使用我提供的函数(玻尔兹曼 sigmoid)强制拟合通过点 (-70,0) 在数学上是不可能的。另一方面,@Cleb 和@BenBolker 解释了如何强制拟合通过任何其他点,例如 (-50, 0.09)。
正如我们在您的问题下方的评论中所讨论的那样,使用您提供的函数(没有偏移量)是不可能强制拟合通过 0 的。
但是,您可以通过为单个数据点设置 weights
来强制曲线通过其他数据点。所以例如如果您给数据点 A 的权重等于 1,数据点 B 的权重等于 1000,则数据点 B 对拟合更重要(就将要最小化的残差总和的贡献而言)比 A 和拟合因此将被迫通过 B。
这是完整的代码,我将在下面进行更详细的解释:
# your data
xdata <- c(-60, -50, -40, -30, -20, -10, -0, 10)
ydata <- c(0.04, 0.09, 0.38, 0.63, 0.79, 1, 0.83, 0.56)
plot(xdata, ydata, ylim=c(0, 1.1))
fit <-nls(ydata ~ 1 / (1 + exp((xdata - p1) / p2)), start=list(p1=mean(xdata), p2=-5))
# plot the fit
xr = data.frame(xdata = seq(min(xdata), max(xdata), len=200))
lines(xr$xdata, predict(fit, newdata=xr))
# set all weights to 1, do the fit again; the plot looks identical to the previous one
we = rep(1, length(xdata))
fit2 = nls(ydata ~ 1 / (1 + exp((xdata - p1) / p2)), weights=we, start=list(p1=mean(xdata) ,p2=-5))
lines(xr$xdata, predict(fit2, newdata=xr), col='blue')
# set weight for the data point -30,0.38, and fit again
we[3] = 1000
fit3 = nls(ydata ~ 1 / (1 + exp((xdata - p1) / p2)), weights=we, start=list(p1=mean(xdata), p2=-5))
lines(xr$xdata, predict(fit3, newdata=xr), col='red')
legend('topleft', c('fit without weights', 'fit with weights 1', 'weighted fit for -40,0.38'),
lty=c(1, 1, 1),
lwd=c(2.5, 2.5, 2.5),
col=c('black', 'blue', 'red'))
输出如下;如您所见,拟合现在通过所需的数据点(红线):
所以这是怎么回事:我首先像你一样进行拟合,然后我进行权重拟合,其中所有权重都设置为 1;因此,该图看起来与之前的相同,蓝线隐藏了黑线。然后 - 对于 fit3
- 我将第三个数据点的权重更改为 1000,这意味着现在最小二乘拟合的 "important" 比其他点多得多,并且新拟合会遍历此数据点(红线)。
这也是我更改行的第二个示例
we[3] = 1000
到
we[2] = 1000
强制拟合通过第二个数据点:
如果您想获得有关 weights
参数的更多信息,您可以在此处阅读:documentation
基于@Cleb 的回答,这里有一种方法可以选择函数必须通过的指定点并求解其中一个参数的结果方程:
dd <- data.frame(x=c(-60,-50,-40,-30,-20,-10,-0,10),
y=c(0.04, 0.09, 0.38, 0.63, 0.79, 1, 0.83, 0.56))
初始拟合(为方便起见,使用 plogis()
而不是 1/(1+exp(-...))
):
fit <- nls(y ~ plogis(-(x-p1)/p2),
data=dd,
start=list(p1=mean(dd$x),p2=-5))
现在插入 (x3,y3)
并求解 p2:
y3 = 1/(1+exp((x-p1)/p2))
logit(x) = qlogis(-x) = log(x/(1-x))
e.g. plogis(2)=0.88 -> qlogis(0.88)=2
qlogis(y3) = -(x-p1)/p2
p2 = -(x3-p1)/qlogis(y3)
设置函数并插入 p2
:
p2 <- function(p1,x,y) {
-(x-p1)/qlogis(y)
}
fit2 <- nls(y ~ plogis(-(x-p1)/p2(p1,dd$x[3],dd$y[3])),
data=dd,
start=list(p1=mean(dd$x)))
绘制结果:
plot(y~x,data=dd,ylim=c(0,1.1))
xr <- data.frame(x = seq(min(dd$x),max(dd$x),len=200))
lines(xr$x,predict(fit,newdata=xr))
lines(xr$x,predict(fit2,newdata=xr),col=2)
我正在尝试将 Boltzmann sigmoid 1/(1+exp((x-p1)/p2))
拟合到这个小型实验数据集:
xdata <- c(-60,-50,-40,-30,-20,-10,-0,10)
ydata <- c(0.04, 0.09, 0.38, 0.63, 0.79, 1, 0.83, 0.56)
我知道这很简单。例如,使用 nls
:
fit <-nls(ydata ~ 1/(1+exp((xdata-p1)/p2)),start=list(p1=mean(xdata),p2=-5))
我得到以下结果:
Formula: ydata ~ 1/(1 + exp((xdata - p1)/p2))
Parameters:
Estimate Std. Error t value Pr(>|t|)
p1 -33.671 4.755 -7.081 0.000398 ***
p2 -10.336 4.312 -2.397 0.053490 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1904 on 6 degrees of freedom
Number of iterations to convergence: 13
Achieved convergence tolerance: 7.079e-06
然而,我需要(由于理论原因)拟合曲线恰好通过点(-70, 0)
。虽然上面显示的拟合表达式的值在 x = -70
处接近零,但它并不完全为零,这不是我想要的。
所以,问题是:有没有办法告诉 nls
(或其他一些函数)拟合相同的表达式,但强制它通过指定的点?
更新:
正如评论中提到的那样,使用我提供的函数(玻尔兹曼 sigmoid)强制拟合通过点 (-70,0) 在数学上是不可能的。另一方面,@Cleb 和@BenBolker 解释了如何强制拟合通过任何其他点,例如 (-50, 0.09)。
正如我们在您的问题下方的评论中所讨论的那样,使用您提供的函数(没有偏移量)是不可能强制拟合通过 0 的。
但是,您可以通过为单个数据点设置 weights
来强制曲线通过其他数据点。所以例如如果您给数据点 A 的权重等于 1,数据点 B 的权重等于 1000,则数据点 B 对拟合更重要(就将要最小化的残差总和的贡献而言)比 A 和拟合因此将被迫通过 B。
这是完整的代码,我将在下面进行更详细的解释:
# your data
xdata <- c(-60, -50, -40, -30, -20, -10, -0, 10)
ydata <- c(0.04, 0.09, 0.38, 0.63, 0.79, 1, 0.83, 0.56)
plot(xdata, ydata, ylim=c(0, 1.1))
fit <-nls(ydata ~ 1 / (1 + exp((xdata - p1) / p2)), start=list(p1=mean(xdata), p2=-5))
# plot the fit
xr = data.frame(xdata = seq(min(xdata), max(xdata), len=200))
lines(xr$xdata, predict(fit, newdata=xr))
# set all weights to 1, do the fit again; the plot looks identical to the previous one
we = rep(1, length(xdata))
fit2 = nls(ydata ~ 1 / (1 + exp((xdata - p1) / p2)), weights=we, start=list(p1=mean(xdata) ,p2=-5))
lines(xr$xdata, predict(fit2, newdata=xr), col='blue')
# set weight for the data point -30,0.38, and fit again
we[3] = 1000
fit3 = nls(ydata ~ 1 / (1 + exp((xdata - p1) / p2)), weights=we, start=list(p1=mean(xdata), p2=-5))
lines(xr$xdata, predict(fit3, newdata=xr), col='red')
legend('topleft', c('fit without weights', 'fit with weights 1', 'weighted fit for -40,0.38'),
lty=c(1, 1, 1),
lwd=c(2.5, 2.5, 2.5),
col=c('black', 'blue', 'red'))
输出如下;如您所见,拟合现在通过所需的数据点(红线):
所以这是怎么回事:我首先像你一样进行拟合,然后我进行权重拟合,其中所有权重都设置为 1;因此,该图看起来与之前的相同,蓝线隐藏了黑线。然后 - 对于 fit3
- 我将第三个数据点的权重更改为 1000,这意味着现在最小二乘拟合的 "important" 比其他点多得多,并且新拟合会遍历此数据点(红线)。
这也是我更改行的第二个示例
we[3] = 1000
到
we[2] = 1000
强制拟合通过第二个数据点:
如果您想获得有关 weights
参数的更多信息,您可以在此处阅读:documentation
基于@Cleb 的回答,这里有一种方法可以选择函数必须通过的指定点并求解其中一个参数的结果方程:
dd <- data.frame(x=c(-60,-50,-40,-30,-20,-10,-0,10),
y=c(0.04, 0.09, 0.38, 0.63, 0.79, 1, 0.83, 0.56))
初始拟合(为方便起见,使用 plogis()
而不是 1/(1+exp(-...))
):
fit <- nls(y ~ plogis(-(x-p1)/p2),
data=dd,
start=list(p1=mean(dd$x),p2=-5))
现在插入 (x3,y3)
并求解 p2:
y3 = 1/(1+exp((x-p1)/p2))
logit(x) = qlogis(-x) = log(x/(1-x))
e.g. plogis(2)=0.88 -> qlogis(0.88)=2
qlogis(y3) = -(x-p1)/p2
p2 = -(x3-p1)/qlogis(y3)
设置函数并插入 p2
:
p2 <- function(p1,x,y) {
-(x-p1)/qlogis(y)
}
fit2 <- nls(y ~ plogis(-(x-p1)/p2(p1,dd$x[3],dd$y[3])),
data=dd,
start=list(p1=mean(dd$x)))
绘制结果:
plot(y~x,data=dd,ylim=c(0,1.1))
xr <- data.frame(x = seq(min(dd$x),max(dd$x),len=200))
lines(xr$x,predict(fit,newdata=xr))
lines(xr$x,predict(fit2,newdata=xr),col=2)