ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

プログラムによる近似

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

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

多項式近似関数

関数

説明

polyfit

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

polyval

polyval(p,x) は、polyfit で決まる n 次多項式の x での値を計算します。

この例では、多項式を用いてデータをモデル化する方法を示します。

時間 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];
plot(t,y,'o')
title('Plot of y Versus t')

2 次の多項式関数を用いて、このデータのモデル化を試みることができます。

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

未知の係数 $a_{0}$, $a_{1}$, and $a_{2}$ は、モデルからデータの偏差の二乗和を最小化することで計算されます (最小二乗近似)。

polyfit を使用して多項式係数を求めます。

p = polyfit(t,y,2)
p =
   -0.2942    1.0231    0.4981

MATLAB は、多項式係数を降べきの順で計算します。

データについての 2 次多項式モデルは、以下の式で与えられます。

$$y = -0.2942t^{2}+1.0231t+0.4981.$$

等間隔の時間 t2 で多項式を計算します。次に、元のデータとモデルを同じ場所にプロットします。

t2 = 0:0.1:2.8;
y2 = polyval(p,t2);
figure
plot(t,y,'o',t2,y2)
title('Plot of Data (Points) and Model (Line)')

データの時間ベクトルでモデルを評価します。

y2 = polyval(p,t);

残差を計算します。

res = y - y2;

残差をプロットします。

figure, plot(t,res,'+')
title('Plot of the Residuals')

2 次の近似はデータの基本的な形に概ね従っていますが、データが乗っているように見える滑らかな曲線をとらえていないように見えることに注意してください。残差にはパターンがあるように見え、他のモデルが必要である可能性を示しています。(次に示す) 5 次多項式は、データが変動するに従って、より適切な近似をします。

今度は polyfit の 5 次多項式を用いて演習を繰り返します。

p5 = polyfit(t,y,5)
p5 =
    0.7303   -3.5892    5.4281   -2.5175    0.5910    0.6000

t2 で多項式を計算し、新規の Figure ウィンドウでデータ上に近似をプロットします。

y3 = polyval(p5,t2);
figure
plot(t,y,'o',t2,y3)
title('Fifth-Degree Polynomial Fit')

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

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

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

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

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

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

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

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
a =

    1.3983
   -0.8860
    0.3085

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

$$y = 1.3983 - 0.8860e^{-t} + 0.3085te^{-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')

重回帰

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

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

$x_{1}$ $x_{2}$ の複数の値で、量 $y$ を測定します。これらの値を、ベクトル x1x2y にそれぞれ格納します。

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
a =

    0.1018
    0.4844
   -0.2847

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

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

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

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

    0.0038

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

プログラムによる近似

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

人口調査データのサンプルを census.mat から読み込みます。このファイルには、1790 年から 1990 年までの米国の人口データが含まれています。

load census

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

  • cdate は、1790 年から 1990 年までの 10 年ごとの年度を示す列ベクトルです。

  • pop は、cdate の中の各年度に対応する米国の人口数が記された列ベクトルです。

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

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

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

相関係数の計算

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

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

corrcoef(cdate,pop)
ans =

    1.0000    0.9597
    0.9597    1.0000

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

データへの多項式近似

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

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

[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)');

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

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

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

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

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

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

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

次のコードでは、±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

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

この情報は役に立ちましたか?