2016-02-03 99 views
1

我想在Matlab中將一些數據擬合成形式爲y = r^n /(r^n + K^n)的Hill函數。我有r,y的數據,我需要找到K,n。無法擬合非線性曲線數據在Matlab中

我在廣泛閱讀文檔後嘗試了兩種不同的方法 - 一種使用曲線擬合工具箱中的fit,另一種使用優化工具箱中的lsqcurvefit。我還沒有取得任何成功。我錯過了什麼?

這是我的外部數據和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]; 

'飛度' 代碼:

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

錯誤消息:

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' 代碼:

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 using snls (line 47) 
Objective function is returning undefined values at initial point. lsqcurvefit cannot 
continue. 
+1

開始通過更換'K.^N'與'K2'。你以後總是可以拿第n根來恢復原來的常量,但是簡化表達式以適應融合。 –

+0

大部分問題可能來自於你使用'x^a'這樣的表達式,並帶有負數'x'(參見'xdata'的前半部分)和部分'a'。這可能是你爲什麼會引用「imag」的錯誤。但除此之外:爲什麼傳遞一個包含函數句柄而不是句柄本身的字符串呢?從[fit]文檔(http://www.mathworks.com/help/curvefit/fit.html#inputarg_fitType):「爲了適應自定義模型,使用MATLAB表達式,線性模型術語的單元陣列,一個匿名函數,或者用'fittype'函數創建'fittype'並將其用作'fitType'參數。「 –

+1

[你可以避免在匿名函數中使用參數數組],它應該是非常簡單的。然而,你的負指數模型並不重要。 –

回答

3

首先,我認爲你需要3個變量,從開始的,因爲山上的功能將是1 max和你的數據最大爲10.所以要麼通過做ydata=ydata./max(ydata)來規範你的數據,要麼添加第三個變量(我只是爲了演示)。這是我做的:

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) 

旁註:

你的情況,你真正想要做的是通過看Y=1./ydata,而不是給你轉換數據,然後進行擬合,然後再採用另一個1./Y以獲得前面「山」函數表示形式的答案。這是因爲你的問題在本質上是雙線性,因此通過去1./ydata你雙線性關係,因而它的順序1 polyfit會做:

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

enter image description here