Solved – 4th parameter of Boltzmann sigmoid must be greater than .9 in R

curve fittingnlsr

I'm trying to fit a 4 parameter boltzmann sigmoid and get an error: "Error in nls(y ~ a0 + (a1 – a0)/(1 + exp((a2 – x)/a3)), start = list(a0 = max(y), :
singular gradient"

I have figured out that the code runs if the a3 parameter is set to .9 or greater, or -.9 or less. Does anyone have the reason this is? I want to provide a starting parameter for a3 as the slope according to the description on this website: http://www.originlab.com/doc/Origin-Help/Boltzmann-FitFunc . That is why I have the linear fit coefficient a3.s, but the result is < .9 and I get the error. Is there a way to estimate a3.s prior to use as starting parameter for nls? I am simply using a linear fit of the midpoint of the sigmoid +/- 10 x units – is that the correct interpretation of the a3 parameter?

Here is my code:

#fit boltzman sigmoid
a0.s=max(y); a1.s=min(y); a2.i=which.min( abs(((a0.s+a1.s)/2) - y) ); a2.s=x[a2.i]
lin.x.i=x<a2.s+10 & x>a2.s-10
a3.s=unname(coef(lm(y[lin.x.i]~so[lin.x.i]))[2])
fit <- nls(y ~ a0 + (a1-a0)/(1+exp((a2-x)/a3)), 
start=list(a0=max(y), a1=min(y), a2=a2.s,a3=.9) , trace=TRUE)
params=coef(fit)
curve(params[1]+(params[2]-params[1])/(1+exp((params[3]-x)/params[4])), 1,100,col='black',add=T,type='l')

Here is the data:

x=c( 75,  40,  90,  55, 15, 100,  10,  70,  90,  50,  15,   5,   5,  70, 100,  20,  60,  65,  20,  50,  30,  85,  60,  80,  55,  40,  45,  95,  10,  55, 60,  10,  35,  80,  75,  25,  30,   5,  35,  50, 100,  40,  30,  80,  20,  45,  25,  25,  95,  95,  65,  35,  90,  85,  70,  15,  75,  45,  85,  65);

y=c(4.673686, 0.034781, 5.014355, 0.843847, 0.013337, 4.214557, 0.015299, 5.017280, 4.327815, 0.041139, 0.008704, 0.007437, 0.005125, 4.725786, 3.869776, 0.018725, 4.514051, 3.232932, 0.012979, 0.257651, 0.028170, 4.723512, 2.676991, 5.018232, 0.633399, 0.040133, 0.051864, 5.019395, 0.006505, 0.642376, 2.752317, 0.010827, 0.029303, 4.050711, 3.698887, 0.018385, 0.029491, 0.013894, 0.032034, 0.053761, 5.029349, 0.038272, 0.032619, 5.030450, 0.022356, 0.053421, 0.025370, 0.024763, 4.948973, 3.254528, 1.149153, 0.038530, 4.612227, 4.048692, 4.809153, 0.016246, 5.014711, 0.062841, 5.026961, 2.951881)

Related to this question: the formula on the linked to website has a slight variation in the equation, where the the a2 parameter is used in the form "x-a2", while the equation I provided, and got from my data acquisition software's curve fitting function is the one I provided in the code with "a2-x". Which form of the Boltzmann is correct? Does the difference matter?

Best Answer

Why are you reinventing the wheel? Use the native function SSfpl for your model.

fitr <- nls(y ~ SSfpl(x, a1, a0, a2, ma3))

The parameter ma3 is -a3 in your notation, but otherwise the parametrization is identical, and you get slightly better convergence.

You should probably be using weighted least squares, since ordinary least squares assumes the variability of $y-E(y)$ does not depend on the value of $x$; which is clearly violated in your data.

Related Solutions

Solved – Using nls() function in R for exponential function

Non-linear least squares solves $min_\beta \sum (y_i-f(x_i;\beta))^2$. This is quadratic in $\beta$ if $f$ is linear in $\beta$. Your $f$ is not linear in $\beta$, so the NLS objective function is not quadratic in $\beta$. Of course, you don't need the function to be quadratic to guarantee convergence to a unique minimum, rather you need $min_\beta \sum (y_i-f(x_i;\beta))^2$ to be convex in $\beta$. Presumably, with your $f$, the NLS objective function is not convex. It doesn't look, to me, like the kind of $f$ which generates a convex objective function. That's pretty much the explanation. You can have lots of minima or one minimum.

If I were fitting the function that you are, I would use an entirely different approach. I would not just blindly use NLS. If you look carefully at your function, $f(x_i;\beta)=a*x_i^2exp(-bx_i)+c$ it is almost linear in the parameters. If you fixed $b$ at some value, say 0.1, then you could fit $a$ and $c$ by OLS: \begin{align} y_i &= a*x_i^2exp(-0.1x_i)+c \\ &= a*z_i+c \end{align} The variable $z_i$ is defined $z_i=x_i^2exp(-0.1x_i)$. This means that, once you have picked $b$, the optimal value of $a=\widehat{Cov}(y,z)/\hat{V}(z)$ and the optimal value of $c=\overline{y}-a*\overline{z}$.

So what, right? At the very least, this is how you should pick starting values for $a$ and $c$. But, really, this reduces the search for optimal parameters to a one dimensional search over $b$. With a modern computer, one dimensional searches are fast and easy. If you have some idea of what reasonable values for $b$ are, then you can just define an interval $[b_{low},b_{high}]$ and grid search for the b which gives the lowest sum of squared errors. Then use that $b$ and its associated optimal $a$ and $c$ to start NLS from.

Or, you could do something more sophisticated. Suppose you are searching over $b$, using the optimal $a(b)$ and $c(b)$ from OLS. Then the NLS objective function is $\sum \left(y_i - f(x_i;a(b),b,c(b))\right)^2$. The envelope theorem makes the derivative of this very easy to calculate: \begin{align} \frac{d}{d b} \sum \left(y_i - f(x_i;\beta)\right)^2 &= \sum 2\left(y_i - f(x_i;\beta)\right)\frac{d}{d b}f(x_i;\beta)\\ &= \sum 2\left(y_i - f(x_i;\beta)\right)(-abx_i^2exp(-bx_i)) \end{align}

So, you can easily write a function to calculate the NLS objective function for any given $b$ and you can easily write a function to calculate the derivative of the NLS objective function for any $b$. These two ingredients are enough to get a optimizer going on your function. Then, after you find the optimal $b$, just run NLS with that $b$ and its associated optimal $a$ and $c$. It will converge in one iteration.

Solved – Logistic growth curve with R nls

Have you tried using SSLogis in your nls call? Right now, you're just fitting a line, and the reason you're getting that error is because nls requires a symbolic variable in the passed formula.


data <- data.frame(
  x = 0:15,
  y = c(3.493, 5.282, 6.357, 9.201, 11.224, 12.964, 16.226, 18.137,
        19.590, 21.955, 22.862, 23.869, 24.243, 24.344, 24.919, 25.108)
)


model = nls(y ~ SSlogis(x, a, b, c), data = data)

plot(data$x, data$y)
lines(data$x, predict(model))

This results in the following fit

Best Answer

Related Solutions

Solved – Using nls() function in R for exponential function

Solved – Logistic growth curve with R nls

Related Question