曲線擬和函數lsqcurvefit&nlinfit

lsqcurvefit(非線性最小二乘法)
help了一下, 發現官方的文檔過於詳盡, 節選部分吧.
LSQCURVEFIT solves non-linear least squares problems.

X = LSQCURVEFIT(FUN,X0,XDATA,YDATA) starts at X0 and finds coefficients X to best fit the nonlinear functions in FUN to the data YDATA (in the least-squares sense).

FUN accepts inputs X and XDATA and returns a vector (or matrix) of function values F, where F is the same size as YDATA, evaluated at X and XDATA.

NOTE: FUN should return FUN(X,XDATA) and not the sum-of-squares sum((FUN(X,XDATA)-YDATA).^2).

((FUN(X,XDATA)-YDATA) is squared and summed implicitly in the algorithm.)

這里面對函數的各個參數的意義和用法都做了解釋.
[x,resnorm]=lsqcurvefit(fun,x0,xdata,ydata,...)
fun 是我們需要擬合的函數，這是重點
x0 是我們對函數中各參數的預測值，這也是重點
xdata 則是橫軸坐標的值
ydata 是縱軸的值
需要注意的是lsqcurvefit里面的fun是我們需要擬和的函數, 需要另外編寫.

下面用以前某年的CUMCM真題修改作為示例.

--------------------------------------------------------------------------------------------
tj/s | 100 | 200 | 300 | 400 | 500 | 600 | 700 | 800 | 900 | 1000 |
--------------------------------------------------------------------------------------------
Cj/(mg*cm-3) | 4.54 | 4.99 | 5.35 | 5.65 | 5.90 | 6.10 | 6.26 | 6.39 | 6.50 | 6.59 |
--------------------------------------------------------------------------------------------
在已知部分參數的情況下, 求函數

求函數的最下值的點(K,a,b)

1. 編寫M文件(curvefun.m)

function f = curvefun(x,tdata)
f = x(1)+x(2)*exp(-0.02*x(3)*tdata);
end

2. 編寫程序(test1.m)

tdata = linspace(100,1000,10);
cdata = 1e-05.*[454 499 535 565 590 610 626 639 650 659];
x0 = [0.2 0.05 0.05];
x = lsqcurvefit(@curvefun, x0, tdata, cdata);
f = curvefun(x, tdata)
plot(tdata, cdata, '*', tdata, f, 'r-')

3. 輸出結果
x= 0.0069 -0.0029 0.0809
即表示k=0.2542, a=0.0063, b=-0.0034

nlinfit

從matlab給出的幫助文檔來看, nlinfit與lsqcurvefit同屬與非線性最小二乘擬和, 一般來說都是能得到比較接近的結果.
但是由於nlinfit使用的是牛頓方法, 在使用是需要給出你和參數的假設初值, 有些問題對初值比較敏感, 不同的初值會導致差異比較大.
下面示例nlinfit的用法:

混凝土的抗壓強度隨着養護時間的延長而增加, 現將一批混凝土作成12個試塊, 記錄了養護日期x(日)以及抗壓強度y(kg/cm2)的數據

-------------------------------------------------------------------------------------------------------
養護時間x | 2 | 3 | 4 | 5 | 7 | 9 | 12 | 14 | 17 | 21 | 28 | 56 |
-------------------------------------------------------------------------------------------------------
抗壓強度y | 35+r | 42+r | 47+r | 53+r | 59+r | 65+r | 68+r | 73+r | 76+r | 82+r | 86+r | 99+r |
--------------------------------------------------------------------------------------------------------
建立非線性回歸模型, 對得到的模型和系數進行檢驗.
注明：此題中的+r代表加上一個[-0.5,0.5]之間的隨機數

1. 編寫m文件(myfunc.m)

function f = myfunc(beta, x)
f = beta(1)+beta(2)*exp(beta(4)*x)+beta(3)*exp(-beta(4)*x);
end

2. 編寫程序

clc;clear;
x=[2 3 4 5 7 9 12 14 17 21 28 56];
r=rand(1,12)-0.5; %rand產生的是0,1之間的隨機數, 這里表示產生12個[-0.5 0.5]之間的隨機數
y1=[35 42 47 53 59 65 68 73 76 82 86 99];
y=y1+r
beta=nlinfit(x,y,@myfunc,[0.5 0.5 0.5 0.5]); %nlinfit的語法與lsqcurvefit基本類似, 只是參數的順序上有些差異, 這里不再贅述
a=beta(1),k1=beta(2),k2=beta(3),m=beta(4)
%test the model
xx=min(x):max(x);
yy=a+k1*exp(m*xx)+k2*exp(-m*xx);
plot(x,y,'o',xx,yy,'r')

當然, 這里如果不想另外編寫m函數, 可以使用上一篇筆記中的內聯函數inline來實現. 只需要將程序里面的beta=.....替換成下面兩行即可, 注意替換后的nlinfit不需要使用@符號了.

myfunc=inline('beta(1)+beta(2)*exp(beta(4)*x)+beta(3)*exp(-beta(4)*x)','beta','x');
beta=nlinfit(x,y,myfunc,[0.5 0.5 0.5 0.5]);

3. 輸出結果

a = 88.0591, k1 = 0.0318, k2 = -62.7924, m = 0.1044

曲線擬和函數lsqcurvefit&nlinfit

免責聲明！