How to fit a nonlinear function with parameter-dependent constraint?

Hello,

I want to fit the nonlinear function

    y(x)=a(1)*(C1/x)^a(2)

to experimental data. Here, a(1) and a(2) are the parameters to be optimized. C1 is a known constant.

In order to avoid singularity at x=0, I'd like to manipulate the fitting function and want to set

    y(x=0)=a(1)*(a(2)+1)*(C1/C2)^a(2)

where C2 ist another known constant. Since y(x=0) depends on a(1) and a(2), I don't see any way to implement this.

Do you have any suggestions?

Thanks in advance!

Answer 1

[Alan Weiss]

Why not take a logarithm?

 log y(x) = log(a(1)) + a(2)*log(C1/x)

Have log(a(1)) and a(2) be the variables to find, you now have a linear regression. Solve it using LinearModel.fit or just plain \ (mldivide).

In detail, your data points are:

 Independent variable: measurements of log(C1/x)
 Response variable:    measurements of log(y)

Make a matrix XX with two columns, the independent variable in column 2, and ones in column 1. Make a matrix YY with the response variable.

Your model is XX*a = YY.

Least squares solution: ahat = XX \ YY.

Alan Weiss

MATLAB mathematical toolbox documentation

[Markus]

Thanks for you post, Alan. Taking the logarithm is a good idea, I obviously forgot about.

However, I still don't see how to implement the function manipulation at x=0 that I stated in the question. Do I overlook something in your suggestion?

[Alan Weiss]

Do you have data with x = 0? If so, is y = Inf there? If x = 0 and y is not = Inf, then your model is no good, and you need to come up with a better model. If you have data with x = 0 and y = Inf, then simply throw those points away, they add nothing to the fitting of the parameters.

If you have no data with x = 0, then don't worry about it.

Alan Weiss

MATLAB mathematical toolbox documentation

[Markus]

To answer your questions: I do have data at x=0 which is not equal to Inf. And yes, my model is not matching the envelope of my experimental data very well, but still its the one I need to fit.

Inspired by your post I solved the problem with "fmincon". I added a third coloumn to XX, with its first row element equal to one, all others equal to zero. Additionally, I set

    p(3)=log10(a(2)+1);

so my parameters are:

    p(1)=log10(a(1)); 
    p(2)=a(2);
    p(3)=log10(a(2)+1);

My modelfunction is:

    f_temp=XX*p;
    f=sqrt(sum(log10(YY)-f_temp)^2);

where f ist my objective function I aim to minimize with fmincon like

    x=fmincon(@modelfunction,x0,[],[],[],[],[],[],@coneq,opt);

with a nonlinear constraint to p(3) given by

    function [c,ceq]=coneq(x)
      c=[];
      ceq=log10(1+p(2))-p(3);
    end

This seems to work for me.