第二十八章

第二十八章

曲線擬合與迴歸分析

本 章 重 點

曲線擬合（Curve Fitting）是資料分析的一個重要

步驟 ， 其 目 的 是 要 經 由 有 限 的 取 樣 點（ Sample

Points）來建立一個數學模型，並藉由此模型來進

行進一步的預測與分析。本章將介紹曲線擬合及迴

歸分析的基本方法，並介紹如何以 MATLAB 來進

行實作。

28-1 線性迴歸

通常在曲線擬合的問題中，所建立的數學模型是「單輸入／單輸出」

（Single-input Single-output，簡稱SISO），所以其特性可用一條曲線 來表示。若欲建立具有兩個輸入的數學模型，則其特性可用一個曲面 來表示，此類問題可稱為曲面擬合（Surface Fitting）。相同的概念，

可以延伸到一般的「多輸入／單輸出」（MISO） 或「多輸入／多輸 出」（MIMO）的數學模型，例如類神經網路（Artificial Neural Networks）

等。

無論是曲線擬合或是曲面擬合（或是其它多輸入模型的擬合問題），

在資料分析上都稱為迴歸分析（Regression Analysis），或稱為資料擬 合（Data Fitting），其牽涉到的數學理論與分析技巧相當廣泛。但在 本章中，我們只介紹基本的方法，以及如何使用 MATLAB 來實現這 些基本的方法。

迴歸分析與所使用的數學模型有很大的關係，如果所使用的模型是線 性模型，則此類問題稱為線性迴歸（Linear Regression）；反之，若使 用非線性模型，則稱為非線性迴歸（Nonlinear Regression）。本節將 介線性迴歸，下節則介紹非線性迴線。

線性迴歸是在迴歸分析中最常用的方法，其解法已建立良久，且其相 關數學性質在一般教科書中都找得到。

假設我們的觀察資料是美國自 1790 至 1990 年（以 10 年為一單位）的 總人口，此資料可由載入檔案 census.mat 得到，如下：

>>load census.mat % 載入人口資料

>> plot(cdate, pop, 'o'); % cdate 代表年度，pop 代表人口總數

17500 1800 1850 1900 1950 2000 50

100 150 200 250

其中 cdate 為年度，pop 為人口總數，兩者都是長度為 21 的向量。假 設我們要預測在西元 2000 年及 2010 年的美國人口總數，我們就必須 根據上述資料來建立一數學模型，並依此模型來進行預測。

由上圖可觀察得知，通過這些點的曲線可能是一條二次的拋物線，所 以我們可以假設所需的數學模型為

2 2 1 0 2 1 0

, , )

;

( x a a a a a x a x f

y = = + +

其中 y（人口總數）為此模型的輸出，x（年度）為此模型的輸入，

a

₀^、

a

1 ^及

a

₂ 則為此模型的參數（Parameters）。由於這些參數相對於輸出 y 是呈線性關係，所以此 模型稱為「具有線性參數 （ Linear-in-the- parameters）」的模型，我們的任務則是找出最好的參數值，使得模型 輸 出 與 實 際 資 料 越 接 近 越 好 ， 此 過 程 即 稱 為 線 性 迴 歸 （ Linear Regression）。

假設我們的觀察資料可寫成

( x

_i

, y

_i

)

, i = 1 ~ 21。當輸入為

x

_i ^時，實際

輸出為

y

_i^{，但模型的預測值為}

f ( x

_i

; a

₀

, a

₁

, a

₂

) = a

₀

+ a

₁

x

_i

+ a

₂

x

_i^2，因此

平方誤差為

[ y

i

− f ⁽ x

i

⁾ ]

^{2，而總平方誤差為：}

[ ] _∑ [ ( ) ]

∑

= =

+ +

−

=

−

=

²¹

1

2 2 2 1 0 21

1

)

2

(

i

i i i

i

i

i

f x y a a x a x

y E

上述平方誤差 E 是參數

a

₀^、

a

₁^、

a

₂ 的函數，因此我們可以求出

E

^對

a

0^、

a

₁^、

a

₂ 的導式，令其為零，再解出

a

₀^、

a

₁^、

a

₂^{。由於此模型具}

線性參數，所以平方誤差

E

^為

a

₀^、

a

₁^、

a

₂ ^{的二次式，而導式}

a

0

E

∂

∂

、

a

1

E

∂

∂

及

a

2

E

∂

∂

為

a

₀^、

a

₁^、

a

₂ 的一次式，因此在令導式為零之後，可以

解出參數

a

₀^、

a

₁ ^及

a

₂ ^{的最佳值。}

從另一方面來看，假設 21 個觀察點均通過此拋物線，則：

21 2 21 2 21 1 0

2 2 2 2 2 1 0

1 2 1 2 1 1 0

y x a x a a

y x a x a a

y x a x a a

= +

+

= +

+

= +

+

M

亦可寫成

{ 1 2 3 M 4

4 3 4

4 2 1

M

A y

y y y

x x

x x

x x

 

 





 

 





=

 







 







 

 







 

 







21 2 1

3 2 1

2 21 21

2 2 2

2 1 1

1 1

1

θ

θ θ θ

其中

A

^、

y

^為已知，

θ

為未知向量。很明顯的，上述方程組含 21 個 方程式，但卻只有 3 個未知數（

θ = [ θ

1

, θ

2

, θ

3

]

^{T），所以解通常不存在，}

我們只能找到一組

θ

，使得等號兩邊的差異為最小，此差異可寫成

) (

)

2

(

θ θ

θ y A y A

A

y − = −

^T

−

此即為前述的總平方誤差

E

^。

由於此類問題經常碰到，所以 MATLAB 提供一簡單方便的「左除」

（\）來解出最佳的

θ

^，如下：

>> A = [ones(size(cdate)), cdate, cdate.^2];

>> y = pop;

>> theta = A\y;

>> theta(1)

ans =

2.1130e+004

>> theta(2)

ans =

-23.5097

>> theta(3)

ans =

0.0065

提 示 ：

左除的概念，可記憶如下：

8

原先的方程式是 A*theta = y，我們可將 A 移項至等號右邊，而得到 theta = A\y。必須小心的是：原先 A 在乘式的第一項，所以移到等號右邊後，A 仍然 必須是除式的第一項。

8

若我們要解的方程式是 theta*A = y，則同樣的概念可得到最小平方解 theta = A/ y。

8

有關「左除」的詳細說明，可見本書第二十四章「線性代數」的第四節「聯立 線性方程式」。

一旦找出最佳的

θ

^{，即可繪圖如下：}

>> plot(cdate, pop, 'o', cdate, A*theta, '-');

>> legend('Actual value', 'Predicted value');

>> xlabel('Year');

>> ylabel('Population');

17500 1800 1850 1900 1950 2000

50 100 150 200 250

Year

Population

A ctual value P redicted value

由此可見到所找到的模型

2 2

2 1

0

21130 23 . 51 0 . 00654 )

( x a a x a x x x

f

y = = + + = − +

已經能夠逼近所給的 21 個觀察點。根據此模型，我們可預測美國在 2000 年的人口總數為：

>> [1,2000,2000^2]*theta % 在 2000 年美國人口線數預測值

ans =

274.6221

而在 2010 年的預測人口總數為：

>> [1,2010,2010^2]*theta % 在 2010 年美國人口線數預測值

ans =

301.8240

在上述例子中，我們的模型是一個二次拋物線，我們可將之推廣，得 到一個 n 次多項式：

n n

x a x

a a x f

y = ( ) =

₀

+

₁

+ L +

利用多項式擬合來進行曲線擬合，通稱為「多項式擬合（Polynomial Fitting）」，也是在資料分析常用到的技巧。針對多項式擬合，MATLAB 提供了 polyfit 指令，其效果和「左除」（\）是一致的，如下例：

>> theta = polyfit(cdate, pop, 2); % 進行多項式擬合，cdate 代 % 表年度，pop 代表人口總數

>> theta(1)

ans =

0.0065

>> theta(2)

ans =

-23.5097

>> theta(3)

ans =

2.1130e+004

其中 2 代表的了用到的模型是 2 次多項式。有關 polyfit 及相關的多項 式，可見本書第二十五章「多項式的處理及分析」。

線性迴歸的成功與否，與所選取的模型有很大的關係。模型所含的參 數越多，平方誤差會越小，若參數個數等於資料點個數，則在一般情 況下，平方誤差會等於零，但這並不表示預測會最準，因為資料點含 有雜訊，完全吻合資料的模型亦代表此模型受雜訊的影響最大，預測 之準確度也會較差。因此，「模型複雜度」（即可變參數的個數）和

「預測準確度」是相互抗衡的兩個因素，無法兩者兼得，因此我們必 須以對問題本身的瞭解來找到一個平衡點。

一般「多輸入／單輸出」的線性迴歸數學模型可寫成

γ χ χ

χ

χ ) ( ) ( ) ( )

( a

₁

f

₁

a

₂

f

₂

a

_n

f

_n

f

y = = + + L +

其中

χ

為輸入（長度為 m 的向量），y 為輸出（純量），

a

₁^、

a

₂^{、… 、}

a

n 為可變的未知參數，

f

i

( χ )

、i = 1~n，則是已知的函數，稱為基底 函數（Basis Functions）。假設所給的資料點為

( χ

_i

, y

_i

), i = 1 L m

^，這

些資料點稱為取樣資料（Sample Data）或訓練資料（Training Data），

將這些資料點帶入模型後可得：

) ( )

( )

( )

(

) ( )

( )

( )

(

2 2 1

1

2 2 1

1 1

m m

m m

1 1

1 1

χ χ

χ χ

χ χ

χ χ

n n m

n n

f a f

a f

a f

y

f a f

a f

a f

y

+ + +

=

=

+ + +

=

=

L M

L

或可表示成矩陣格式：

{ {

y m n

A

m n m

n

y y a

a f

f

f f

 







 







=

 







 







 







 







M M

4 4 4

4 3

4 4 4

4 2

1 L

M O M

L

_{1 1}

1

1 1

1

) ( )

(

) ( )

(

θ

x x

x x

由於在一般情況下，m >> n（即資料點個數遠大於可變參數個數），

因此上式無解，因此欲使上式成立，須加上一誤差向量 e：

y e è + = A

平方誤差則可寫成

) (

) (

)

( θ e

²

e e y A θ y A θ

E = =

^T

= −

^T

−

而使

E ( θ )

^為最小的

θ

值可用 MATLAB 的「左除」來算出，即

y

= A\

θ

^。

提 示 ：

8

理論上，最佳的 θ ^值為 ( ) ^A

^T

^A

⁻¹

^A

^T

^{y ，但是} (

^ATA

)

⁻¹

^{容易造成電腦內部計}

算的誤差，MATLAB 的「左除」可以得到較穩定且正確的數值解。

8

欲知如何推導最 θ 值，可參考筆者另一著作：「Neuro – Fuzzy and Soft Computing」，Prentice Hall，1997。

以下以“ peaks ”函數為例，來說明一般的線性迴歸。若在 MATLAB 命此函數可表示成：

2 2 2

2 2

2 ( 1) 3 5 ( 1)

2

3 1 10 5

) 1 (

3

^{x y}

x x y e

^{x y}

e

^{x y}

e x

z

^{− − +}



^{− −}

−

^{− + −}



 



 − −

−

−

=

在此例中，我們假設：

l 基底函數已知

l 訓練資料包含正規分佈的雜訊

因此上述函數可寫成：

n y x f y x f y x f

n e

y x f y x f

n e

e y x x

e x z

y x

y x y

x y

x

+ +

+

=

+

−

−

=

+

 −



 



 − −

−

−

=

− +

−

− +

−

−

− +

−

−

) , ( )

, ( )

,

( 3

) 1 , ( 10 ) , ( 3

3 1 10 5

) 1 ( 3

3 3 2

2 1

1

) 1 ( 2

1

) 1 ( 5

3 )

1 ( 2

2 2

2 2 2

2 2

2

θ θ

θ

其中我們假設

θ

₁^、

θ

2 ^和

θ

3 ^{是未知參數，}

n

則是平均為零、變異為 1 的 正規分佈雜訊。例如：如果要取得 100 筆訓練資料，可輸入：

>> point_n = 10;

>> [xx, yy, zz] = peaks(point_n);

>> zz = zz + randn(size(zz)); % 加入雜訊

>> surf(xx, yy, zz);

>> axis tight

-3 -2 -1 0 1 2 3 -2

0 2 -4 -2 0 2 4 6

在上例中，randn 指令的使用即在加入正規分佈雜訊。

上圖為我們收集到的訓練資料，由於雜訊很大，所以和原先未帶雜訊 的圖形差異很大。現在我們要用已知的基底函數，來找出最佳的

θ

₁^、

θ

2 ^和

θ

3^，如下：

>> x = xx(:);

>> y = yy(:);

>> z = zz(:);

>> A = [(1-x).^2.*exp(-(x.^2)-(y+1).^2), (x/5-x.^3- y.^5).*exp(-x.^2-y.^2), exp(-(x+1).^2-y.^2)];

>> theta = A\z

theta =

3.0794

-9.1175

0.4095

由此找出的

θ

^{值和最佳值}





 



 − −

3 , 1 10 ,

3

相當接近。根據此參數，我

們可以輸入較密的點，得到迴歸後的曲面：

>> point_n = 31;

>> [xx, yy] = meshgrid(linspace(-3, 3, point_n), linspace(-3, 3, point_n));

>> x = xx(:);

>> y = yy(:);

>> A = [(1-x).^2.*exp(-(x.^2)-(y+1).^2), (x/5-x.^3- y.^5).*exp(-x.^2-y.^2), exp(-(x+1).^2-y.^2)];

>> zz = reshape(A*theta, point_n, point_n);

>> surf(xx, yy, zz);

>> axis tight

在上圖中，可知迴歸後的曲面和原先的曲面相當接近。最主要的原因 是：我們猜對了基底函數（或是更正確的說，我們偷看了正確的基底

函數），因此得到非常好的曲面擬合。一般而言，若不知正確的基底 函數而胡亂選用，很難由 3 個可變函數達到 100 個資料點的良好擬 合。

在上例中我們曾在資料點加入正規分佈（Normal Distributed）的雜訊。

事實上，只要基底函數正確，而且雜訊是正規分佈，那麼當資料點越 來越多，上述的最小平方法就可以逼近參數的真正數值。欲詳知此理 論的細節，可詳閱筆者的著作「Neural-Fuzzy and Soft Computing」，

Prentice Hall ，1997。

28-2 ^{非線性迴歸}

相對而言，非線性迴歸（Nonlinear Regression）是一個比較困難的問 題，原因如下：

l 無法一次找到最佳解。

l 無法保證能夠找到最佳解。

l 須引用各種非線性最佳化的方法。

l 各種相關數學性質並不明顯。

以數學來描述，假設所用的數學模型是

( r , θ r ) x f

y =

^，其中

x ^r

^是輸入向

量，

θ r

是可變非線性函數，y 是輸出變數，則總平方誤差為

[ ]

²

1

) , ( )

( ∑

=

−

=

^m

i

i

i

f x

y

E θ r r θ r

其中

( x , r

i

y

i

)

是第 I 個已知資料點。由於

θ r

是

f

^{的非線性參數，所以}

)

( θ r

E

^並不是

θ r

的二次式，因此我們並無法由

( θ r ) E

^對

θ r

的導式為零

來 解 出 最 佳 的

θ r

值 。 退 而 求 其 次 ， 我 們 必 須 用 一 般 最 佳 化

（Optimization）的方法，來找出

( θ r )

E

的最小值，例如梯度下降法

（Gradient Descent），或是 Simplex 下坡式搜尋（Simplex Downhill search）等。

舉例來說，假設所用的數學模型為

x

x

a e

e a

y =

₁ ^λ1

+

₂ ^λ2

其中，

a

₁ ^、

a

₂^{為線性參數，但 λ}1、λ2 為非線性參數，則此模型為 非線性，總平方誤差可表示如下：

( )

∑

=

+

−

=

^m

i

x x

i

a e a e

y a

a

E

ⁱ

1

2 2 1

2 1 2 1

2 2

)

1

, , ,

( λ λ

^{λ λ}

欲找出使

E

^為最小的

a

₁ ^、

a

₂^{、λ1及 λ2}，我們需將 E 寫成一函式，

並 由 其 它 最 佳 化 的 方 法 來 求 出 此 函 式 的 最 小 值 。 假 設 此 函 式 為 fun_e1，其內容可顯示如下：

>> type fun_e1

function E = fun_e1(theta, data) x = data(:,1);

y = data(:,2);

model_y = theta(1)*exp(theta(3)*x)+theta(2)*exp(theta(4)*x);

E = sum((y-model_y).^2);

其中 theta 是參數向量，包含了

a

₁ ^、

a

₂^、λ1及 λ2，data 則是觀察到 的資料點，傳回的值

E

則是總平方誤差。欲求出

E

^{的最小值， 我}

們可使用 fminsearch 指令，例如：

>>data = [0 5.8955;

0.1000 3.5639;

0.2000 2.5173;

0.3000 1.9790;

0.4000 1.8990;

0.5000 1.3938;

0.6000 1.1359;

0.7000 1.0096;

0.8000 1.0343;

0.9000 0.8435;

1.0000 0.6856;

1.1000 0.6100;

1.2000 0.5392;

1.3000 0.3946;

1.4000 0.3903;

1.5000 0.5474;

1.6000 0.3459;

1.7000 0.1370;

1.8000 0.2211;

1.9000 0.1704;

2.0000 0.2636];

>> init_theta = [0 0 0 0];

>> theta = fminsearch('fun_e1', init_theta, [], data);

>> x = data(:, 1);

>> y = data(:, 2);

>> estimated_y =

theta(1)*exp(theta(3)*x)+theta(2)*exp(theta(4)*x);

>> plot(x, y, 'ro', x, estimated_y, 'b-');

>> legend('Sample data', 'Regression curve');

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0

1 2 3 4 5 6

S am ple data R egression curve

上圖的曲線即為 fminsearch 指令所產生的迴歸曲線。在上述程式中，

data 矩陣包含自變數（即data(:,1)）和因變數（即 data(:,2)），以方便 將 之 傳 入 函 式 fun_e1 。 init_theta 則 是 可 變 參 數 theta 的 起 始 值 。 fminsearch 指 令 則 是 一 個 使 用 Simplex 下 坡 式 搜 尋 法 （ Downhill Simplex Search）的最佳化方法，用來找出 fun_e1 的極小值，並傳回 theta 的最佳值。

上述方法是一個一般性的方法，因為它把對所有參數都是一視同仁，

全部視為非線性參數。但是我們還可以進一步改良，也就是將線性與 非線性參數分開，各用不同的方法來處理。以上例而言，數學模型為

x

x

a e

e a

y =

₁ ^λ1

+

₂ ^λ2

其中

a

₁ ^及

a

₂^{為線性參數，λ1及 λ2}為非線性參數，因此我們可以對 它們分開處理，如下：

l 線性參數（

a

₁ ^及

a

₂）：最小平方法，即「左除」或「\」。

l 非線性參數（λ1及 λ2）：Simplex 下坡式搜尋（即 fminsearch）。

此方法的好處是：最小平方法能夠在非線性參數（λ1及 λ2）固定的 情況下，一次找到最好的線性參數（

a

₁ ^及

a

₂^{）的值，因為搜尋空間}

的維度由 4（

a

₁ ^、

a

₂^{、λ1、λ2），降為 2（λ1及 λ2），最佳化會}

更有效率。

如果要使用上述混成（Hybrid）的方法，函式 fun_e1 須改寫成 fun_e2，

如下：

>> type fun_e2

function E = fun_e2(lambda, data) x = data(:,1);

y = data(:,2);

A = [exp(lambda(1)*x) exp(lambda(2)*x)];

a = A\y;

model_y = a(1)*exp(lambda(1)*x)+a(2)*exp(lambda(2)*x);

E=sum((y-model_y).^2);

function E = fun_e2(lambda, data) x = data(:,1);

y = data(:,2);

A = [exp(lambda(1)*x) exp(lambda(2)*x)];

a = A\y;

model_y = a(1)*exp(lambda(1)*x)+a(2)*exp(lambda(2)*x);

E=sum((y-model_y).^2);

function E = fun_e2(lambda, data) x = data(:,1);

y = data(:,2);

A = [exp(lambda(1)*x) exp(lambda(2)*x)];

a = A\y;

model_y = a(1)*exp(lambda(1)*x)+a(2)*exp(lambda(2)*x);

E = sum((y-model_y).^2);

其中 lambda 是非線性參數向量，只包含 λ1 及 λ2，data 仍是觀察到 的資料點，a 則是利用最小平方法算出的最佳線性參數向量（即

a

₁ ^及

a

2^{，隨傳入的 λ1及 λ2而變），而傳回的}

E

^{仍是總平方誤差。此函}

數和　fun_e1　的最大不同在於我們只須傳入線性參數（λ1及 λ2），

最佳線性參數（

a

₁ ^及

a

₂）的值則在函式中由最小平方法決定。

欲用此混成法求出 E 的最小值，可輸入如下：

>> init_lambda = [0 0];

>> lambda = fminsearch('fun_e2', init_lambda, [], data);

Warning: Rank deficient, rank = 1 tol = 2.1368e-014.

> In d:\mlbook\examples\FUN_E2.M at line 5

In C:\MATLABR11\toolbox\matlab\funfun\fminsearch.m at line 103

Warning: Rank deficient, rank = 1 tol = 2.1374e-014.

> In d:\mlbook\examples\FUN_E2.M at line 5

In C:\MATLABR11\toolbox\matlab\funfun\fminsearch.m at line 160

Warning: Rank deficient, rank = 1 tol = 2.1370e-014.

> In d:\mlbook\examples\FUN_E2.M at line 5

In C:\MATLABR11\toolbox\matlab\funfun\fminsearch.m at line 201

Warning: Rank deficient, rank = 1 tol = 2.1371e-014.

> In d:\mlbook\examples\FUN_E2.M at line 5

In C:\MATLABR11\toolbox\matlab\funfun\fminsearch.m at line 216

Optimization terminated successfully:

the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-004

and F(X) satisfies the convergence criteria using OPTIONS.TolFun of 1.000000e-004

>> x = data(:, 1);

>> y = data(:, 2);

>> A = [exp(lambda(1)*x) exp(lambda(2)*x)];

>> a = A\y;

>> estimated_y = a(1)*exp(lambda(1)*x)+a(2)*exp(lambda(2)*x);

>> plot(x, y, 'ro', x, estimated_y, 'b-');

>> legend('Sample data', 'Regression curve');

0 0.5 1 1.5 2

0 1 2 3 4 5 6

S am ple data R egression curve

此混成法顯然較快而且較準。

提 示 ：

8

在上例中，MATLAB 會產生警告訊息，這是由於在最佳化的過程中，可能出現 A 的 rank 小於其行數，導致 a = A\y 出現 「Rank deficient 」的警告訊息。

但只要最後我們檢視擬合圖形，能得到滿意結果，即可忽略這最佳化過程所產 生的暫時警告訊息。

8

欲檢視上述迴歸過程的動畫顯示，可在 MATLAB 指令視窗下輸入 fitdemo。

除此之外，我們亦可利用變形法（Transformation），將一數學模型轉 換成只包含線性參數的模型。例如，假設一模型為：

ae

bx

y =

取自然對數，可得

bx a y = ln + ln

此時，

ln a

^及

b

變成線性參數，我們可用“最小平方法”找出其值。

請特別注意，此最小平方法所得最小化的總平方誤差是

( )

²

1

ln ln ' ∑

=

−

−

=

^m

i

i

i

a bx

y E

而不是原模型的總平方誤差：

( )

∑

=

−

=

^m

i

bx i

ae

i

y E

1

2

通常

E '

^{為最小值時，}

E

不一定是最小值，但亦離最小值不遠矣！若 要精益求精，可再用 fminsearch 求取

E

的最小值，並以最小平方法得 到的 a 及 b 為搜尋的起點。

使用變形法的例子如下：

>> a = 1; b = -2; % 未知的參數

>> x = (0:0.1:1)'; % 已知資料點的 x 座標

>> y = a*exp(b*x)+randn(size(x))/20; % 已知資料點的 y 座標（含雜訊）

>> A = [ones(size(x)) x];

>> b = log(y);

>> theta = A\b;

>> a1 = exp(theta(1)); % 辨識得到之參數

>> b1 = theta(2); % 辨識得到之參數

>> plot(x, y, 'o', x, a1*exp(b1*x));

>> legend('Actual value', 'Predicted value');

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

A ctual value P redicted value