Non-linear data fitting to a model

I am trying to fit some data using the following code:

xdata = [0.0600    0.1250    0.2540    0.4050    0.5450];
ydata = [0.0800    0.1000    0.1200    0.1400    0.1600];
  objfcn = @(x)(x(1)*x(2))./((x(2)*(1-xdata)./xdata)+(1-xdata)*(x(2)-1))-ydata;
  x0 = rand(2,1);
  opts = optimoptions(@lsqnonlin,'Algorithm','levenberg-marquardt','Display','iter');
  [x,resnorm]=lsqnonlin(objfcn,x0,[],[],opts);
  M = x(1)
  B = x(2)

I am getting results but sometimes M value is coming negative which is not desirable. Also, the predicted ydata deviates largely from the original ydata values. If someone could be help to analyze the small code mentioned above and what is the wrong with results.

Thanks,

Answer 1

[Matt J]

Random guesses for the initial point x0 are a bad idea. I generate a better initial guess below, by re-arranging the model equation as linear. Also, if it is undesirable for M to be negative, then you should use the lb argument to bound it. This means you will then have to use the trust-region algorithm instead of Levenberg-Marquardt, since the latter doesn't support bound constraints.

   xdata=xdata(:); ydata=ydata(:);
   objfcn = @(x)(x(1)*x(2))./((x(2)*(1-xdata)./xdata)+(1-xdata)*(x(2)-1))-ydata;

    expr1=(1-xdata)./xdata;
    expr2=(1-xdata);
    A=[ones(size(xdata)),   -(expr1+expr2).*ydata, ];
    b=ydata.*expr2;
    s=A\b;
    res=A*s-b

    x0 = [s(1)/s(2),s(2)],
    opts = optimoptions(@lsqnonlin,'Display','iter');
    [xf,resnorm]=lsqnonlin(objfcn,x0,[0,-inf],[],opts);

This leads to

      K>> yfit= objfcn(xf)+ydata

    yfit =

        0.0827
        0.0903
        0.1067
        0.1342
        0.1757

I tend to think this is the best fit that will be possible given such a small number of xdata, ydata samples.

[KB]

Thank you Matt, this helps. Now, just curious to know can someone make a reasonable initial guess to have better results ? Actually, later I tried the similar thing using lower bound (trust-region) as you said but with random guess and a different initial guess, however the results were never close to the experimental value of ydata. As you have made a nice initial guess, the outcome is very good. I am just wondering is there any proper way to do make initial guess ? What are the xf value are you getting ? I am trying to run the similar code what you posted but it says initial value is undefined. s(1)/s(2) is coming 0/0.

[Matt J]

Generating the initial guess is a problem-specific art and not a science, I'm afraid. A frequent strategy, though, is to ignore measurement noise and then re-arrange or transform the model equation to something simpler. When we ignore noise in your original model, it can be re-arranged into the equation

 (x(1)*x(2))=((x(2)*(1-xdata)./xdata)+(1-xdata)*(x(2)-1)).*ydata

Making the change of variables s(1)=x(1)*x(2), s(2)=x(2) this becomes linear in s

s(1)-((s(2)*(1-xdata)./xdata)+(1-xdata)*(s(2)-1)).*ydata = 0

and can be solved using standard linear solvers, which is what I did for you in my code above.

Another way to generate an initial guess, because you only have 2 unknown variables, is to make a surf plot of the resnorm function and inspect it for the approximate location of the minimum. Or maybe instead of plotting, just sample the resnorm function on some coarse sampling grid and use min() to locate an approximate global minimum.