最小二乘法曲線擬合

閱讀原文

預備知識 最小二乘法, Nelder-Mead 演算法

   我們在 「最小二乘法」 中見到的三個例子中, 方差函數都是待定係數的線性組合, 這種情況下我們令偏導為零後得到的是線性方程組, 便於求解. 然而當方差不是待定係數的線性組合時, 得到的方程組往往非常複雜, 這時就需要藉助數值計算. 相比用數值計算解 N 元的非線性方程組, 更簡單的方法是直接用數值方法尋找方差函數的極小值(如 Nelder-Mead 演算法) . 實踐證明, 大多數情況下極小值點僅有一個, 那就是最小值點.

   為了驗證結果的正確性, 我們先來用數值方法擬合 Acos(x+φ0)+C , 並與「最小二乘法」 中的方法比較結果.

圖1:運行結果

curveFit.m

0001 function curveFit0002 close all;0003 % 隨機生成簡諧曲線0004 N = 20;0005 x0 = linspace(0, 2*pi, N);0006 y0 = 5*rand * sin(x0 + 2*pi * rand) + 10 * rand;0007 y0 = y0 + 2*rand(1,20); % 產生隨機誤差0008 scatter(x0, y0, +); % 畫出散點0009 hold on;0010 0011 % Nelder-Mead 求方差最小值點0012 f = @(x) norm( x(1)*sin(x0 + x(2)) + x(3) - y0 );0013 c = NelderMead(f, [5, 1, 5], 1e-7);0014 % 畫擬合結果0015 x = linspace(0, 2*pi, 50);0016 y1 = c(1) * sin(x + c(2)) + c(3);0017 plot(x, y1);0018 0019 % 解線性方程組求方差最小值點0020 c = sinfit(x0, y0);0021 % 畫擬合結果0022 y2 = c(1)*cos(x) + c(2)*sin(x) + c(3);0023 plot(x, y2, .);0024 end0025 0026 % 擬合簡諧曲線0027 % y = C(1)*cos(x) + C(2)*sin(x) + C(3)0028 function C = sinfit(x, y)0029 N = numel(x);0030 cosx = cos(x); sinx = sin(x);0031 sc = sum(sinx.*cosx);0032 s = sum(sinx); c = sum(cosx);0033 % 係數矩陣0034 M = [sum(cosx.^2), sc , c;0035 sc, sum(sinx.^2), s;0036 c, s, N];0037 b = [sum(y.*cosx); sum(y.*sinx); sum(y)];0038 C = Mb; % 解線性方程組0039 end

(剩下部分見頂部的「閱讀原文」)

推薦閱讀:

TAG:演算法 | 計算物理學 | 數值分析 |