Search code examples
matlabcurve-fittingnonlinear-functions

Unable to fit nonlinear curve to data in Matlab


I am trying to fit some data in Matlab to a Hill function of the form y = r^n/(r^n+K^n). I have data for r,y and I need to find K,n.

I tried two different approaches after reading the docs extensively - one uses fit from the CurveFitting Toolbox and the other uses lsqcurvefit from the Optimization Toolbox. I haven't had any success with either. What am I missing?

Here is my xdata and ydata:

xdata = logspace(-2,2,101);
ydata = [0.0981 0.1074 0.1177 0.1289 0.1411 0.1545 0.1692 0.1852 0.2027 ...
          0.2219 0.2428 0.2656 0.2905 0.3176 0.3472 0.3795 0.4146 0.4528 ...
          0.4944 0.5395 0.5886 0.6418 0.6994 0.7618 0.8293 0.9022 0.9808 ...
          1.0655 1.1566 1.2544 1.3592 1.4713 1.5909 1.7183 1.8537 1.9972 ...
          2.1490 2.3089 2.4770 2.6532 2.8371 3.0286 3.2272 3.4324 3.6437 ...
          3.8603 4.0815 4.3065 4.5344 4.7642 4.9950 5.2258 5.4556 5.6833 ...
          5.9082 6.1292 6.3457 6.5567 6.7616 6.9599 7.1511 7.3347 7.5105 ...
          7.6783 7.8379 7.9893 8.1324 8.2675 8.3946 8.5139 8.6257 8.7301 ...
          8.8276 8.9184 9.0029 9.0812 9.1539 9.2212 9.2834 9.3408 9.3939 ...
          9.4427 9.4877 9.5291 9.5672 9.6022 9.6343 9.6638 9.6909 9.7157 ...
          9.7384 9.7592 9.7783 9.7957 9.8117 9.8263 9.8397 9.8519 9.8630 ...
          9.8732 9.8826];

'Fit' code:

HillEqn = '@(x,xdata)xdata.^x(1)./(xdata.^x(1)+x(2).^x(1))';
startPoints = [1 1];
fit(xdata',ydata',HillEqn,'Start',startPoints)

Error message:

Error using fittype>iTestCustomModelEvaluation (line 726)
Expression @(x,xdata)xdata.^x(1)./(xdata.^x(1)+x(2).^x(1)) is not a valid MATLAB
expression, has non-scalar coefficients, or cannot be evaluated:
Undefined function 'imag' for input arguments of type 'function_handle'.

'lsqcurvefit' code:

fun = @(x,xdata) xdata.^x(1)./(xdata.^x(1)+x(2).^x(1));
x0 = [1,1]; % Initial Parameter Estimates
x = lsqcurvefit(fun,x0,xdata,ydata);

Error message:

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

Solution

  • First I think you need 3 variables to start from, because the hill function will be max of 1, and your data it is maxed at 10. So either normalize your data by doing ydata=ydata./max(ydata), or add a 3rd variable (which I did just for the demonstration). This is how I did it:

    startPoints = [1 1 1];
    s = fitoptions('Method','NonlinearLeastSquares',... %
        'Lower',[0    0   0  ],...
        'Upper',[inf  inf inf],...
        'Startpoint',startPoints);
    
    HillEqn = fittype( 'x.^a1./(x.^a1+a2.^a1)*a3','options',s);
    [ffun,gofrr] = fit(xdata(:),ydata(:),HillEqn);
    yfit=feval(ffun,xdata(:)); %Fitted function
    plot(xdata,ydata,'-bx',xdata,yfit,'-ro');
    

    enter image description here

    ffun = 
    
     General model:
     ffun(x) = x.^a1./(x.^a1+a2.^a1)*a3
     Coefficients (with 95% confidence bounds):
       a1 =       1.004  (1.004, 1.004)
       a2 =      0.9977  (0.9975, 0.9979)
       a3 =       9.979  (9.978, 9.979)
    

    Side note:

    In your case what you really want to do is to transform you data by looking at Y=1./ydata instead, then fit, and then take another 1./Y to get the answer in the previous "hill" function representation. This is because your problem is bilinear in nature , so by going 1./ydata you get a bilinear relation, for which a polyfit of order 1 will do:

    Y=1./ydata;
    X = 1./xdata;
    p=polyfit(X,Y,1);
    plot(X,Y,'-bx',X,polyval(p,X),'-ro')
    

    enter image description here