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.
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');
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')