non linear regression problem. fitnlm gives error

Hello

I need to do curve fitting with 2 set of data and I know the model function and initial guess. I used "fitnlm " but it gives error.

would you please check my code and tell me whether I made a mistake or it is the nature of these data.

load FITDATAK

load FITDATAB

tb1= table (FITDATA1,FITDATA)

modelfun= @(b,x) b(1)*x(:,1)^b(2)

beta0= [0.22 -0.5]

md1= fitnlm(tb1,modelfun,beta)

it gives me this errors :

Error using nlinfit>checkFunVals

The function you provided as the MODELFUN input has returned Inf or NaN values.

Error in nlinfit>LMfit (line 596)

if funValCheck && ~isfinite(sse), checkFunVals(r); end

Error in nlinfit (line 284)

[beta,J,~,cause,fullr] = LMfit(X,yw, modelw,beta,options,verbose,maxiter);

Error in NonLinearModel/fitter (line 1153)

nlinfit(X,y,F,b0,opts,wtargs{:},errormodelargs{:});

Answer 1

[Torsten]

Look at the values of K and B. Do you really want to approximate B by a*K^b ?

format long
K = load("FITDATAK.mat");
K = K.FITDATA1

K = 526×1

1.0e+00 * 0.000000000165736 0.000000000146554 0.000000000138203 0.000000000104830 0.000000000067225 0.000000000053905 0.000000000043810 0.000000000022345 0.000000000019131 0.000000000009961

nnz(isnan(K))

ans = 
     0

nnz(isinf(K))

ans = 
     0

B = load("FITDATAB.mat");
B = B.FITDATA

B = 526×1

308500 439800 410200 749100 721800 900400 1072400 1304600 3641600 4397100

nnz(isnan(B))

ans = 
     0

nnz(isinf(B))

ans = 
     0

modelfun = @(b) b(1)*K.^b(2)-B;
beta0= [1e17  1];
sol = lsqnonlin(modelfun,beta0,[],[],optimset('MaxFunEvals',100000,'MaxIter',100000))

Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

sol = 1×2

1.0e+17 * 0.999999999999999 0.000000000000000

norm(modelfun(sol))

ans = 
     6.834566371856458e+10

[sepideh]

Dear Torsten

I am sure about the model function. it has been proved analytically that B and K has the relationship :

there is no doubt that K has a negative exponential but there is a huge controversy around the value of c2. there is a range for c2 from (-0.5 to -1.5) or even more, among experimental data.

but dimentional analysis insists that the correct exponential is -0.5 (which I put as an initial guess)

we decided to take a universal fit of the data to support the analytical result (-0.5)

I did it by excel but the resultant fit is not robust and acts strangely:

here is the result of power-law fitting in excel while the axis are in logarithmic scale

then I check an alternative way in which I do the fitting with logarithmic data and linear fit :

results are same but the correlation declined!

I want to check this fit with matlab too to see the difference

[Torsten]

K = load("FITDATAK.mat");
K = K.FITDATA1;
B = load("FITDATAB.mat");
B = B.FITDATA;

% Fit log(B) against log(b(1)) + b(2)*log(K)
M = [ones(numel(K),1) log(K)];
rhs = log(B);
pars = M\rhs;
b1 = exp(pars(1))
b2 = pars(2)
% Plot result
loglog(K,B,'o')
hold on
kmin = min(K);
kmax = max(K);
k = linspace(kmin,kmax,1000);
loglog(k,b1*k.^b2)

% Use result of linearized fit as initial guess for nonlinear fit and solve
modelfun = @(b,xdata)b(1)*xdata.^b(2);
fun = @(b,xdata) modelfun(b,xdata)-B;
sol = lsqnonlin(@(b)fun(b,K),[b1 b2],[],[],optimset('MaxFunEvals',100000,'MaxIter',100000))
% Plot result
loglog(k,modelfun(sol,k))
hold off

[sepideh]

when I use the following code

load logk

load logbeta

x=logk\logbata

MATLAB generate the similar answer but I dont know if I can trust this result or not. so I decide to use a non liner regression

[Torsten]

My code from above gives the result of the linearized and the nonlinear problem in one plot.

The result for the nonlinear problem formulation looks far off, but remember that the plot is in loglog scale where the result of the linearized problem is the optimal one.

[sepideh]

Perfect and precise

thank you

Answer 2

[Matt J]

Your model function can easily generate NaN's and Infs because with fitnlm there is nothing to bound the b parameters (see below).. It would be better to use fit() with the 'power1' model, and with appropriate bounds on b(2).

load FITDATAK

modelfun= @(b,x) b(1)*x^b(2);

modelfun([1,-100],FITDATA1(50))

ans = Inf

modelfun([0,-100],FITDATA1(50))

ans = NaN