Friday, December 18, 2009

Nonlinear least squares help

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';