プログラムによる近似

多項式モデルのための 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. One or more of the lines displays its values using only markers

重回帰

ライブスクリプトを開く

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

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 a line object which displays its values using only markers.

プロットは、変数間に高い相関があることを表す強いパターンを示しています。

相関係数の計算

例では、ここで、統計的な相関を決め、データのモデリングが適切であるか確かめるために、変数 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, xlabel Census Year, ylabel Population (millions) contains 2 objects of type line. One or more of the lines displays its values using only markers 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 a line object which displays its values using only markers.

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

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

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

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

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

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

[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, xlabel Census Year, ylabel Population (millions) contains 4 objects of type line. One or more of the lines displays its values using only markers

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