Sometimes it can help to have different solvers to debug a non-linear least squares problem. The R function nls did not return a solution, but CONOPT clearly showed the regression model was not very good. An improved model is proposed.

Afterwards an even better model was contributed by another poster.

$ontext

*Problem: *

*I was trying to estimate the weibull model using nls after putting OLS*

*values as the initial inputs to NLS.*

*I tried multiple times but still i m getting the same error of Error in*

*nlsModel(formula, mf, start, wts) :*

*singular gradient matrix at initial parameter estimates. *

*The Program is as below *

*> vel <- c(1,2,3,4,5,6,7,8,9,10,11,12,13,14)*

*> df <- data.frame(conc, vel)*

*> df*

*conc vel*

*1 0.077 1*

*2 0.328 2*

*3 0.882 3*

*4 1.195 4*

*5 1.884 5*

*6 3.577 6*

*7 6.549 7*

*8 13.000 8*

*9 33.690 9*

*10 52.220 10*

*11 90.140 11*

*12 166.050 12*

*13 233.620 13*

*14 346.890 14*

*> plot(df$vel, df$conc)*

*> para0.st <- c(K=450,*

*+ alpha=0.054,beta=3.398 )*

*> fit0 <- nls(*

*+ conc~ K-(K*exp(-(vel/alpha)^beta)), df, start= para0.st,trace=T)*

*Error in nlsModel(formula, mf, start, wts) :*

*singular gradient matrix at initial parameter estimates *

*Indeed: this model does not seem not to be suited for the data. WITH CONOPT I can solve the*

*least squares problem but the results are: *

*---- 102 VARIABLE r.L residuals *

*i1 -67.787, i2 -67.536, i3 -66.982, i4 -66.669, i5 -65.980, i6 -64.287, i7 -61.315*

*i8 -54.864, i9 -34.174, i10 -15.644, i11 22.276, i12 98.186, i13 165.756, i14 279.026 *

*---- 102 VARIABLE sse.L = 150223.167 sum of squared errors*

*VARIABLE K.L = 67.864 to be estimated*

*VARIABLE alpha.L = 0.054 to be estimated*

*VARIABLE beta.L = 3.398 to be estimated *

*A better fit results with the model conc = K + alpha*vel^beta *

*---- 125 VARIABLE r.L residuals *

*i1 0.833, i2 1.073, i3 1.541, i4 1.503, i5 1.174, i6 0.469, i7 -1.466, i8 -4.083*

*i9 1.077, i10 -5.453, i11 -6.091, i12 12.766, i13 -1.379, i14 -1.964 *

*---- 125 VARIABLE sse.L = 263.633 sum of squared errors*

*VARIABLE K.L = -0.756 to be estimated*

*VARIABLE alpha.L = 2.816336E-4 to be estimated*

*VARIABLE beta.L = 5.317 to be estimated *

*Regression results for this model are: *

*Parameter Estimate Std. Error t value Pr(>|t|)*

*K -7.55855E-01 1.88999E+00 -3.99926E-01 6.96867E-01*

*alpha 2.81634E-04 1.11730E-04 2.52066E+00 2.84420E-02 **

*beta 5.31695E+00 1.51652E-01 3.50602E+01 1.22025E-12 *** *

*In R we can do now: *

*> vel <- c(1,2,3,4,5,6,7,8,9,10,11,12,13,14)*

*> conc <- c(0.077,0.328,0.882,1.195,1.884,3.577,6.549,13.000,33.690,52.220,90.140,166.050,233.620,346.890)*

*> df<-data.frame(conc,vel)*

*> nls(conc~K+alpha*vel^beta,df,start<-c(K=-0.7,alpha=0.0003,beta=5.3),trace=T)*

*349.1687 : -0.7000 0.0003 5.3000*

*264.4474 : -0.7507993362 0.0002808334 5.3172851705*

*263.6331 : -0.7559514801 0.0002816438 5.3169303333*

*263.6331 : -0.7558506044 0.0002816331 5.3169456595*

*Nonlinear regression model*

*model: conc ~ K + alpha * vel^beta*

*data: df*

*K alpha beta*

*-0.7558506 0.0002816 5.3169457*

*residual sum-of-squares: 263.6 *

*Number of iterations to convergence: 3*

*Achieved convergence tolerance: 1.599e-06 *

*An even better model was suggested: conc = alpha*exp[beta/vel]: *

*Parameter Estimate Std. Error t value Pr(>|t|)*

*alpha 3.57740E+04 4.18811E+03 8.54180E+00 1.91040E-06 ****

*beta -6.50221E+01 1.54097E+00 -4.21955E+01 2.03681E-14 *** *

*with sse: 249.6362 *

*See: http://permalink.gmane.org/gmane.comp.lang.r.general/173730 *

$offtext

**set** i /i1*i14/;

**table** data(i,*)

conc vel

i1 0.077 1

i2 0.328 2

i3 0.882 3

i4 1.195 4

i5 1.884 5

i6 3.577 6

i7 6.549 7

i8 13.000 8

i9 33.690 9

i10 52.220 10

i11 90.140 11

i12 166.050 12

i13 233.620 13

i14 346.890 14

;

**parameters** conc(i),vel(i);

conc(i) = data(i,'conc');

vel(i) = data(i,'vel');

**variable**

sse 'sum of squared errors'

r(i) 'residuals'

K 'to be estimated'

alpha 'to be estimated'

beta 'to be estimated'

;

**equations**

calcsse

nlfit(i) 'non-linear fit'

;

K.L=450;

alpha.L=0.054;

beta.L=3.398;

calcsse.. sse=e=**sum**(i,sqr(r(i)));

nlfit(i).. conc(i) =e= K-(K*exp(-[vel(i)/alpha]**beta)) + r(i);

**model** m /calcsse,nlfit/;

**option** nlp=conopt;

**solve** m minimizing sse using nlp;

**display** r.L,SSE.L,K.L,alpha.L,beta.L;

***

** the next model shows a better fit*

** *

**equation** nlfit2(i) 'non-linear fit';

nlfit2(i).. conc(i) =e= K+alpha*[vel(i)**beta] + r(i);

**model** m2 /calcsse,nlfit2/;

**solve** m2 minimizing sse using nlp;

**display** r.L,SSE.L,K.L,alpha.L,beta.L;

**option** nlp=nls;

**solve** m2 minimizing sse using nlp;

**parameter** r1(i); r1(i) = r.l(i);

**scalars** a1,b1; a1=alpha.l; b1=beta.l;

***

** even better (suggested by http://permalink.gmane.org/gmane.comp.lang.r.general/173730)*

** *

**equation** nlfit3(i) 'non-linear fit';

nlfit3(i).. conc(i) =e= alpha*exp[beta/vel(i)] + r(i);

**model** m3 /calcsse,nlfit3/;

**option** nlp=conopt;

**solve** m3 minimizing sse using nlp;

**display** r.L,SSE.L,K.L,alpha.L,beta.L;

**option** nlp=nls;

**solve** m3 minimizing sse using nlp;

**parameter** r2(i); r2(i) = r.l(i);

**scalars** a2,b2; a2=alpha.l; b2=beta.l;

***

** plot results*

***

**file** pltdat /exp.dat/;

**loop**(i,

**put** pltdat vel(i):17:5,conc(i):17:5,r1(i):17:5,r2(i):17:5/;

);

**putclose**;

**file** plt /exp.plt/;

**put** plt;

**put** "K=",K.l:0:16/;

**put** "a1=",a1:0:16/;

**put** "b1=",b1:0:16/;

**put** "f1(x)=K+a1*x**b1"/;

**put** "a2=",a2:0:16/;

**put** "b2=",b2:0:16/;

**put** "f2(x)=a2*exp(b2/x)"/;

**putclose** 'set term png size 600,800'/

'set output "exp.png"'/

'set key left'/

'set multiplot layout 2,1'/

'plot "exp.dat" using 1:2 title "data",f1(x) title "f1(x)=k+a*x**b",f2(x) title "f2(x)=a*exp(b/x)"'/

'plot "exp.dat" using 1:3 title "residuals f1(x)","exp.dat" using 1:4 title "residuals f2(x)"'/

'unset multiplot'/;

**execute** '=wgnuplot.exe exp.plt';

**execute** '=shellexecute exp.png';