nlinfit with modelfun as an integral

Dear all,

I have discrete data A(x,y), which I want to fit by a specified function y=f(x). My function f(x) has the following restriction: df(x)/dlog(x) is equal to a sum of two Gaussians. df(x)/dlog(x) being the derivative of f(x) with respect to the argument log(x). Given this restriction f(x) can be expressed as an integral. This gives 6 fit parameters, namely: the heights of the Gaussians, their means and their standard deviations. I have gone this far:

par=[1e-9,1e-5,1,10000,1,1]; %initial Gaussian parameter guess
Integrand = @(x) ((par(3)/(par(5)*(2*pi)^0.5))*exp((-(log10(x)-log10(par(1))).^2)/(2*par(5)^2))+(par(4)/(par(6)*(2*pi)^0.5))*exp((-(log10(x)-log10(par(2))).^2)/(2*par(6)^2)))./(x*log(10)); %two gaussians
Integral = @(x) integral(Integrand,0,x);
nlinfit(A(:,1),A(:,2),Integral,par)

I however get the following error:

Error using nlinfit (line 142)
Error evaluating model function '@(x)integral(Integrand,0,x)'.

Caused by:
    Error using @(x)integral(Integrand,0,x)
    Too many input arguments.

How can I fix this? Thank you

Answer 1

[Star Strider]

There are a few problems I can see. There could be more, but this should get you started:

Add ‘par’ to the Integrand function arguments:

Integrand = @(par,x) ((par(3)/(par(5) ...

The Integral function then becomes:

Integral = @(par,x) integral(@(x) Integrand(par,x),0,x);

I suggest that you name your initial parameter estimates ‘par0’ (or something other than ‘par’) to avoid confusion:

par0 = [1e-9,1e-5,1,10000,1,1]; %initial Gaussian parameter guess;
est_par = nlinfit(A(:,1),A(:,2),Integral,par0)

The ‘est_par’ assignment are the parameters estimated by nlinfit.

I can’t test your code (so no guarantees), but these changes should at least allow it to run.

[Armantas]

Thanks!

I have now incorporated your suggestions:

par0=[1e-9,1e-5,1,10000,1,1];

Integrand = @(par0,x) ((par0(3)/(par0(5)*(2*pi)^0.5))*exp((-(log10(x)-log10(par0(1))).^2)/(2*par0(5)^2))+(par0(4)/(par0(6)*(2*pi)^0.5))*exp((-(log10(x)-log10(par0(2))).^2)/(2*par0(6)^2)))./(x*log(10));

Integral = @(par0,x) integral(@(x) Integrand(par0,x),0,x);
est_par = nlinfit(A(:,1),A(:,2),Integral,par0)

However, now I get the following error:

Error using nlinfit (line 142)
Error evaluating model function '@(par0,x)integral(@(x)Integrand(par0,x),0,x)'.

Caused by:
    Error using integral (line 86)
    A and B must be floating point scalars.

Something still not right.

My A(x,y) comes as columns btw. Could that be causing the error somehow? I have tried changing it to rows, same error... Thank you!

[Star Strider]

I forgot that integral doesn’t take vector limits, since I never used it with vector limits. The solution is to use a loop, but that requires a separate function file.

The function file I attached here works with vector-valued x, and contains your code for ‘Integrand’. I couldn’t use ‘Integral’ as a variable name for obvious reasons (it causes recursion errors), so I changed its name in the function.

The Integral function file takes vectors for x and integrates from 0 to each x as I believe you intend. It returns a column of results in y, so it should work with your values in nlinfit. Note that your (x,y) arguments have to conform to your original column-oriented A(:,1),A(:,2) configuration for it to work.

Other than adding the loop to allow it to use vector-valued x, it’s entirely your code. I didn’t change anything else.

[Armantas]

Thank you, it now works but nlinfit returns complex valued parameters. That is unphysical given my model i.e. should not be allowed. I googled a bit but couldn't find an easy solution. Could you please point me in the right direction so I can add this restriction?

[Star Strider]

My pleasure!

The only problem I can see is with log10(x), log10(par(1)) and log10(par(2)). You have to fit only x>0, and if you are using lsqcurvefit, constrain par(1) and par(2) to each be >0 (setting a lower bound of eps or perhaps 100*eps for them is the easiest way to accomplish that).

I didn’t see anything else that could be giving you complex values for either your parameters or your function evaluations. I can’t run your code, so this is based on reading it and evaluating Integral with par0 and a vector of x. With your initial parameter estimates and x>0, it returned only real values when I tested it, although the values only changed for small values of x and appeared to be constant (asymptotic to about 10001) for values of x larger than about 1.

Consider:

x = linspace(1E-10,1E-3);
par0 = [1e-9,1e-5,1,10000,1,1];
Q = Integral(par0, x);

figure(1)
plot(x,Q)
grid

[Armantas]

It now works fine! I simply removed the x<0 values from my data vector A(x,y), which I want to fit.

Yes, it would be comfortable to incorporate the bounds into the fit function itself. I have thus also tried your suggestion with the lsqcurvefit by replacing the nlinfit function but I get an error... My code looks more or less like this now:

par=[1.3E-10 1.5E-8 5.3E-5 9E-5 0.63 1.01];   %these fit my data using nlinfit nicely

fit_par = lsqcurvefit(@(par,x)Integral(par,x),par,x,y,1e-14,1e-4);   1e-14 and 1e-4 being my lower and upper x bounds.

%fit_par = nlinfit(x,y,@(par,x)Integral(par,x),par);   used previously

The error message:

Error using snls (line 48)
Objective function is returning undefined values at initial point. lsqcurvefit cannot continue.

Error in lsqncommon (line 175)
    [xC,FVAL,LAMBDA,JACOB,EXITFLAG,OUTPUT,msgData]=...

Error in lsqcurvefit (line 249)
[xCurrent,Resnorm,FVAL,EXITFLAG,OUTPUT,LAMBDA,JACOB] = ...

How can I avoid it?

[Star Strider]

I don’t know what is throwing the error. ‘Undefined values’ I assume to be NaN or ±Inf, so you have to evaluate your function with your initial values of x and parameters to see what the problem is.

However, if you got it to fit well with nlinfit, I would just stop there. Just go with the parameters nlinfit estimated and be happy with them since they seem to provide good results. You can also get confidence intervals on them with nlparci and confidence intervals on the fitted data with nlpredci.

I suggested lsqcurvefit in case you wanted to constrain par(1) and par(2) to be greater than zero to avoid complex results. Since your model obviously converged nicely on positive numbers without constraints with nlinfit, you don’t need to use lsqcurvefit.

................................

With respect to lsqcurvefit though, you can’t set bounds on x — x is your independent variable. You have to have matching x and y values, so if you have an x that is out-of-bounds, you have to delete it and its corresponding y value from the data you give to lsqcurvefit.

You can constrain your parameters with lsqcurvefit. You have 6 parameters so you have to have vectors of upper and lower bounds that are the same dimensions as your parameter vector. So, for instance, if your parameter vector is a (6x1) column vector, your upper and lower bound vectors each have to be (6x1) column vectors as well. There have to be as many constraints as you have parameters. If you don’t want to constrain a particular parameter, set its bounds (in the respective vectors) as -Inf and +Inf. To set a parameter to be greater than zero, set its lower bound to eps.

[Armantas]

Yes, it turns out NaN values in my A(x,y) where causing the problem.

Thank you for all the info, it was very helpful! Now the script is up and running :)

[Star Strider]

My pleasure!

Yours is the most unusual curve-fitting design (with the integral) that I’ve thus far encountered, so I learned much from it. It is also interesting that The Statistics Toolbox function nlinfit could deal with the NaN values on its own.