Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

I'm trying to fit a line+exponential curve to some data. As a start, I tried to do this on some artificial data. The function is: $$y=a+b\cdot r^{(x-m)}+c\cdot x$$ It is effectively an exponential curve with a linear section, as well as an additional horizontal shift parameter (m). However, when I use R's nls() function I get the dreaded "singular gradient matrix at initial parameter estimates" error, even if I use the same parameters that I used to generate the data in the first place.
I've tried the different algorithms, different starting values and tried to use optim to minimise the residual sum of squares, all to no avail. I've read that a possible reason for this could be an over-parametrisation of the formula, but I don't think it is (is it?)
Does anyone have a suggestion for this problem? Or is this just an awkward model?

A short example:

#parameters used to generate the data
reala=-3
realb=5
realc=0.5
realr=0.7
realm=1
x=1:11 #x values - I have 11 timepoint data
#linear+exponential function
y=reala + realb*realr^(x-realm) + realc*x
#add a bit of noise to avoid zero-residual data
jitter_y = jitter(y,amount=0.2)
testdat=data.frame(x,jitter_y)

#try the regression with similar starting values to the the real parameters
linexp=nls(jitter_y~a+b*r^(x-m)+c*x, data=testdat, start=list(a=-3, b=5, c=0.5, r=0.7, m=1), trace=T)

Thanks!

share|improve this question
2  
Hint: look at the coefficient of $r^x$ (for a fixed $r$) and note that $b r^{-m} = \text{constant}$ has a one-dimensional family of solutions $(b,m)$ with $b = r^m \cdot \text{constant}$. –  whuber Jul 14 '11 at 17:35
1  
This isn't an identified model, unless $b$ or $r$ are somehow constrained. I think requiring the $r \in (0,1)$ would do the job. –  Macro Jul 14 '11 at 19:02

2 Answers 2

up vote 3 down vote accepted

I've got bitten by this recently. My intentions were the same, make some artificial model and test it. The main reason is the one given by @whuber and @marco. Such model is not identified. To see that, remember that NLS minimizes the function:

$$\sum_{i=1}^n(y_i-a-br^{x_i-m}-cx_i)^2$$

Say it is minimized by the set of parameters $(a,b,m,r,c)$. It is not hard to see that the the set of parameters $(a,br^{-m},0,r,c)$ will give the same value of the function to be minimized. Hence the model is not identified, i.e. there is no unique solution.

It is also not hard to see why the gradient is singular. Denote

$$f(a,b,r,m,c,x)=a+br^{x-m}+cx$$

Then

$$\frac{\partial f}{\partial b}=r^{x-m}$$

$$\frac{\partial f}{\partial m}=-br^{x-m}$$

and we get that for all $x$

$$b\frac{\partial f}{\partial b}+\frac{\partial f}{\partial m}=0.$$

Hence the matrix

\begin{align} \begin{pmatrix} \nabla f(x_1)\\ \vdots\\ \nabla f(x_n) \end{pmatrix} \end{align}

will not be of full rank and this is why nls will give the the singular gradient message.

I've spent over a week looking for bugs in my code elsewhere till I noticed that the main bug was in the model :)

share|improve this answer

The answers above are, of course, correct. For what its worth, in addition to the explanations given, if you are trying this on an artificial data set, according to the nls help page found at: http://stat.ethz.ch/R-manual/R-patched/library/stats/html/nls.html

R's nls wont be able to handle it. The help page specifically states:

Warning

Do not use nls on artificial "zero-residual" data.

The nls function uses a relative-offset convergence criterion that compares the numerical imprecision at the current parameter estimates to the residual sum-of-squares. This performs well on data of the form

y = f(x, θ) + eps

(with var(eps) > 0). It fails to indicate convergence on data of the form

y = f(x, θ)

because the criterion amounts to comparing two components of the round-off error. If you wish to test nls on artificial data please add a noise component, as shown in the example below.

So, no noise==no good for R's nls.

share|improve this answer
    
Welcome to the site, @B_D_Dubbya. I took the liberty of formatting your answer, I hope you don't mind. You can find more info about editing your answers on CV here. –  gung Sep 13 '12 at 1:54

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.