プログラムによる近似

多項式モデルのための MATLAB 関数

2 つの MATLAB^® 関数が、多項式モデルを用いてデータをモデリングできます。

多項式近似関数

関数	説明
`polyfit`	`polyfit(x,y,n)` は、モデルからのデータの偏差の二乗和を最小にすることにより (最小二乗法)、データ `y` を近似する`n` 次多項式 `p(x)` の係数を求めます。
`polyval`	`polyval(p,x)` は、`polyfit` で決まる `n` 次多項式の `x` での値を計算します。

物理的な状況のモデリングを試みる場合、現在の状況において特定次数のモデルが有意義であるかについて考察することは常に重要です。

非多項式項をもつ線形モデル

ライブスクリプトを開く

この例では、非多項式項を含む線形モデルによって、データを近似する方法を示します。

多項式関数がデータに対して十分なモデルを与えない場合、非多項式項をもつ線形モデルを使用して試すことができます。たとえば、パラメーター $a_{0}$ 、 $a_{1}$ 、 $a_{2}$ に関しては線形で、 $t$ のデータに関しては非線形であるような、次の関数を考えます。

$y = a_{0} + a_{1} e^{- t} + a_{2} t e^{- t} .$

連立方程式を作成して解き、パラメーターを求めることによって、未知の係数 $a_{0}$ 、 $a_{1}$ 、 $a_{2}$ を計算できます。次の構文では、"計画行列" を作成することでこれを行います。ここで、各列は応答 (モデルの項) の予測に使用する変数を表し、各行はそれらの変数の 1 つの観測に対応します。

t と y を列ベクトルとして入力します。

t = [0 0.3 0.8 1.1 1.6 2.3]';
y = [0.6 0.67 1.01 1.35 1.47 1.25]';

計画行列を作成します。

X = [ones(size(t))  exp(-t)  t.*exp(-t)];

モデルの係数を計算します。

a = X\y

したがって、データのモデルは次のように与えられます。

$y = 1.3983 - 0.8860 e^{- t} + 0.3085 t e^{- t} .$

次に、等間隔の点においてモデルを計算し、元のデータとモデルをプロットします。

T = (0:0.1:2.5)';
Y = [ones(size(T))  exp(-T)  T.*exp(-T)]*a;
plot(T,Y,'-',t,y,'o'), grid on
title('Plot of Model and Original Data')

Figure contains an axes object. The axes object with title Plot of Model and Original Data contains 2 objects of type line.

重回帰

ライブスクリプトを開く

この例では、複数の予測子変数が含まれる関数のモデルデータに重回帰を使用する方法を示します。

y が複数の予測子変数を含む関数の場合、変数の関係を表す行列方程式を、データの追加に合わせて拡張しなければなりません。これは、"重回帰" と呼ばれます。

$x_{1}$ と $x_{2}$ のいくつかの値に対して量 $y$ を測定します。これらの値を、ベクトル x1、x2、y にそれぞれ格納します。

x1 = [.2 .5 .6 .8 1.0 1.1]';
x2 = [.1 .3 .4 .9 1.1 1.4]';
y  = [.17 .26 .28 .23 .27 .24]';

データの多変量モデルを、次のように考えます。

$y = a_{0} + a_{1} x_{1} + a_{2} x_{2} .$

重回帰によって、未知の係数 $a_{0}$ 、 $a_{1}$ 、 $a_{2}$ は、モデルからのデータの偏差の二乗和を最小化することで求められます (最小二乗近似)。

計画行列 X を形成することによって、連立方程式を作成して解きます。

X = [ones(size(x1))  x1  x2];

バックスラッシュ演算子を使用してパラメーターを求めます。

a = X\y

データの最小二乗近似モデルは、次のようになります。

$y = 0.1018 + 0.4844 x_{1} - 0.2847 x_{2} .$

モデルを検証するために、モデルから求めた点と観測データ点との偏差の絶対値から最大値を求めます。

Y = X*a;
MaxErr = max(abs(Y - y))

MaxErr = 0.0038

この値がデータ値に比べて十分小さいので、データへの近似がうまくいっているモデルであることを確信できます。

プログラムによる近似

この例では、MATLAB 関数を用いて以下の操作を行う方法を示します。

相関係数の計算
データへの多項式近似
信頼限界のプロットと計算

1790 年から 1990 年までの米国の人口データを含む census.mat から、サンプルの人口調査データを読み込みます。

load census

これにより MATLAB ワークスペースに次の 2 つの変数が追加されます。

cdate は 1790 年から 1990 年までの 10 年ごとの年度を含む列ベクトルです。
pop は cdate の各年に対応する米国の人口数をもつ列ベクトルです。

データをプロットします。

plot(cdate,pop,'ro')
title('U.S. Population from 1790 to 1990')

Figure contains an axes object. The axes object with title U.S. Population from 1790 to 1990 contains an object of type line.

このプロットは、変数間の高い相関を表す明確なパターンを示しています。

相関係数の計算

例では、ここで、統計的な相関を決め、データのモデリングが適切であるか確かめるために、変数 cdate と pop を決定します。相関係数の詳細は、線形相関を参照してください。

相関係数行列を計算します。

corrcoef(cdate,pop)

ans = 2×2

    1.0000    0.9597
    0.9597    1.0000

行列の対角要素は、1 に等しく、各変数がそれ自身と完全に相関していることを表します。非対角要素は 1 に非常に近く、変数 cdate と変数 pop 間に統計上の強い相関があることを示します。

データへの多項式近似

ここでは、MATLAB 関数 polyfit と polyval を適用してデータをモデル化します。

近似パラメーターを計算します。

[p,ErrorEst] = polyfit(cdate,pop,2);

近似を評価します。

pop_fit = polyval(p,cdate,ErrorEst);

データと近似をプロットします。

plot(cdate,pop_fit,'-',cdate,pop,'+');
title('U.S. Population from 1790 to 1990')
legend('Polynomial Model','Data','Location','NorthWest');
xlabel('Census Year');
ylabel('Population (millions)');

Figure contains an axes object. The axes object with title U.S. Population from 1790 to 1990 contains 2 objects of type line. These objects represent Polynomial Model, Data.

このプロットは、2 次の多項式の近似がデータを良く近似していることを示しています。

この近似の残差を計算します。

res = pop - pop_fit;
figure, plot(cdate,res,'+')
title('Residuals for the Quadratic Polynomial Model')

Figure contains an axes object. The axes object with title Residuals for the Quadratic Polynomial Model contains an object of type line.

残差のプロットがパターンを示していることに注意してください。これは、2 次多項式がこのデータのモデル化に適切でない可能性があることを意味します。

信頼限界のプロットと計算

信頼限界は、予測された応答に対する信頼区間です。区間の幅は、近似がどの程度正確であるかを示します。

この例では、サンプルデータ census に polyfit と polyval を適用し、2 次多項式モデルの信頼限界を求めます。

次のコードでは、 $\pm 2 Δ$ の区間を使用します。これは、大規模な標本に対する 95% 信頼区間に相当します。

近似と予測誤差推定 (デルタ) を評価します。

[pop_fit,delta] = polyval(p,cdate,ErrorEst);

データ、近似、信頼限界をプロットします。

plot(cdate,pop,'+',...
     cdate,pop_fit,'g-',...
     cdate,pop_fit+2*delta,'r:',...
     cdate,pop_fit-2*delta,'r:'); 
xlabel('Census Year');
ylabel('Population (millions)');
title('Quadratic Polynomial Fit with Confidence Bounds')
grid on

Figure contains an axes object. The axes object with title Quadratic Polynomial Fit with Confidence Bounds contains 4 objects of type line.

95% の区間は、新しい観測が区間内に入る可能性が 95% であることを示します。