fitting implicit functions with NLINFIT

Is there any way to use implicit functions when fitting them to data with NLINFIT? Details: I need to fit a model to some data with NLINFIT. The problem is that the model to fit is an implicit function, like this:

y= x*(1-exp(a*(h*y-T)))

('x' is the independent variable, 'a' and 'h' are the parameters to be calculated and 'T' is a constant). How can I write a function that returns 'y' when it depends on the values of 'y' itself? I saw a (could-be) solution involving FSOLVE, but I don't have the corresponding toolbox.

Just in case, this comes from integrating (for 't' between 0 and T) the ODE:

dN/dt= a*N /(1+ a*h*N)

given that:

y= N(0)-N(T)

Thanks for any suggestions

Answer 1

[Tom Lane]

NLINFIT wants a response vector Y and a function with unknown parameters. You could try supplying the response vector

Y = zeros(size(y))

and the function

(-y) + x*(1-exp(a*(h*y-T)))

and see if that works. I'd expect it to work. You'd pack x and y together as columns in your X matrix, and unpack them to compute the function value.

[Francisco de Castro]

Thanks Tom, Your suggestion seemed to actually work... however I can't be sure if the solution is 'reasonable' because I get a negative estimate for one of the parameters (which makes no sense). I didn't see any way option to constrain the parameter in NLINFIT. Is there any?

[Star Strider]

The nlinfit function does not allow parameter constraints. Use lsqcurvefit if you need to do that.

Answer 2

[Star Strider]

When I integrate to get N(t) as:

n = dsolve( diff(N) == a*N / (1 + a*h*N), 'IgnoreAnalyticConstraints', true )

in the Symbolic Math Toolbox, the solution is:

N(t) = [0;    lambertw(0, a*h*exp(C3 + a*t))/(a*h)]

and then when I have it create an implicit function for the solution, produces:

Nt =  @(C3,a,h,t)[0.0;  lambertw(0,a.*h.*exp(C3+a.*t))./(a.*h)];

adding N(0) = N0 as an initial condition:

n = dsolve( diff(N) == a*N / (1 + a*h*N), N(0) == N0, 'IgnoreAnalyticConstraints', true )

and integrating yields:

N(t) = lambertw(0, N0*a*h*exp(a*t)*exp(N0*a*h))/(a*h)

and as an anonymous function:

Nt0 = @(N0,a,h,t) lambertw(0,N0.*a.*h.*exp(a.*t).*exp(N0.*a.*h))./(a.*h);

The lambertw function exists in MATLAB (although I admit I've not heard of it until now).

You need to combine parameters N0, a, and h into one vector for curve-fitting purposes in nlinfit or lsqcurvefit. (I do not have the Curve Fitting Toolbox, and use the Statistics and Optimization Toolboxes functions.)

This at least solves your problem of having y on both sides of the equation you want to fit.