Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Welcome To Ask or Share your Answers For Others

Categories

0 votes
994 views
in Technique[技术] by (71.8m points)

r - Forcing nls to fit a curve passing through a specified point

I'm trying to fit a Boltzmann sigmoid 1/(1+exp((x-p1)/p2)) to this small experimental dataset:

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)

I know that it is pretty simple to do it. For example, using nls:

fit <-nls(ydata ~ 1/(1+exp((xdata-p1)/p2)),start=list(p1=mean(xdata),p2=-5))

I get the following results:

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

However, I need (due to theoretical reasons) the fitted curve to pass precisely through the point (-70, 0). Although the value of the fitted expression showed above passes near zero at x = -70, it is not exactly zero, which is not what I want.

So, the question is: Is there a way to tell nls (or some other function) to fit the same expression but forcing it to pass through a specified point?

Update:

As it has been mentioned in the comments, it is mathematically impossible to force the fit to go through the point (-70,0) using the function I provided (the Boltzmann sigmoid). On the other hand, @Cleb and @BenBolker have explained how to force the fit to go through any other point, for instance (-50, 0.09).

See Question&Answers more detail:os

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
Welcome To Ask or Share your Answers For Others

1 Reply

0 votes
by (71.8m points)

It is not possible to force the fit to go through 0 using the function you provide (without an off-set) as we discussed in the comments below your question.

However, you can force the curve to go through other data points by setting weights for individual data points. So e.g. if you give a data point A a weight equals 1 and a data point B a weight equals 1000, the data point B is much more important (in terms of the contribution to the sum of residuals which is going to be minimized) for the fit than A and the fit will therefore be forced to go through B.

Here is the entire code and I explain it in more detail below:

# 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'))

The output looks as follows; as you can see the fit now goes through the desired data point (red line):

enter image description here

So what is going on: I first fit as you did, then I fit with weights whereby all weights are set to 1; therefore, the plot looks identical to the one before and the blue line hides the black line. Then - for fit3 - I change the weight for the third data point to 1000 which means that it is now much more "important" for the least square fit than the other points and the new fit goes through this data point (the red line).

Here is also a second example where I changed the line

we[3] = 1000

to

we[2] = 1000

which forces the fit to go through the second data point:

enter image description here

If you want to get more information about the weights argument you can read here: documentation


与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
OGeek|极客中国-欢迎来到极客的世界,一个免费开放的程序员编程交流平台!开放,进步,分享!让技术改变生活,让极客改变未来! Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Click Here to Ask a Question

...