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';
No comments:
Post a Comment