ドキュメンテーション

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

非線形回帰

パラメトリック非線形回帰モデルとは

パラメトリック非線形モデルは、連続応答変数と複数の連続予測子変数との関係を次の形式で表します。

y = f(X,β) + ε,

ここで

  • y は、応答変数の観測値の n 行 1 列のベクトルです。

  • f は X と β からなる任意の関数で、X の各行と β ベクトルを評価して、y の該当する行の予測値を計算します。

  • X は各行が 1 つの観測値、各列が 1 つの予測子を表す、n 行 p 列の予測子の行列です。

  • β は、推定する p 行 1 列の未知パラメーター ベクトルです。

  • ε は、n 行 1 列の独立な同一の分布に従うランダム外乱ベクトルです。

これに対して、ノンパラメトリック モデルは、モデル パラメーターをもつ予測子と応答の関係を特徴付けようとしません。分類木と回帰木 の場合のように、しばしばグラフィカルな説明になります。

fitnlm は、観測された応答 y とモデル f の予測子 (X,β) の間の平均二乗誤差を最小化する、β パラメーターの値を見つけようとします。そうするには、beta0 を開始値として β ベクトルを最小平均二乗誤差をもつベクトルに修正する手順を反復する必要があります。

データの準備

回帰の近似を開始するには、データを近似関数に望ましい形式にします。すべての回帰手法は、配列 X の入力データと独立したベクトル y の応答データか、テーブルまたはデータセット配列 tbl 内の入力データと tbl の列としての応答データで始まります。入力データの各行が、1 つの観測値を表します。各列が 1 つの予測子 (変数) を表します。

テーブルまたはデータセット配列 tbl では、次のように 'ResponseVar' の名前と値のペアで応答変数を示します。

mdl = fitlm(tbl,'ResponseVar','BloodPressure');

応答変数は既定で最後の列です。

非線形回帰の場合、数値の "カテゴリカル" 予測子は使用できません。カテゴリカル予測子は可能性の固定セットからの値をとります。

入力データと応答データの両方で、欠損値を NaN で表します。

入力と応答データのデータセット配列

たとえば、Excel® スプレッドシートからデータセット配列を作成するには、次のようにします。

ds = dataset('XLSFile','hospital.xls',...
    'ReadObsNames',true);

ワークスペース変数からデータセット配列を作成するには、次のようにします。

load carsmall
ds = dataset(Weight,Model_Year,MPG);

入力および応答データのテーブル

Excel スプレッドシートからテーブルを作成するには、次のようにします。

tbl = readtable('hospital.xls',...
    'ReadRowNames',true);

ワークスペース変数からテーブルを作成するには、次のようにします。

load carsmall
tbl = table(Weight,Model_Year,MPG);

入力データの数値配列と応答の数値ベクトル

たとえば、ワークスペース変数から数値配列を作成するには、次のようにします。

load carsmall
X = [Weight Horsepower Cylinders Model_Year];
y = MPG;

Excel スプレッドシートから数値配列を作成するには、次のようにします。

[X Xnames] = xlsread('hospital.xls');
y = X(:,4); % response y is systolic pressure
X(:,4) = []; % remove y from the X matrix

sex などの数値以外のエントリは X には表示されません。

非線形モデルの表現

非線形モデルを表現するには、いくつか方法があります。最も便利な方法を使用してください。

非線形モデルは、modelfunfitnlm の必須入力となります。

fitnlm は、応答関数 f(X,β) がパラメーター β について滑らかであると仮定します。関数が滑らかでない場合、fitnlm は最適なパラメーター推定値を算出できないことがあります。

無名関数または関数ファイルの関数ハンドル

関数ハンドル @modelfun(b,x) は、ベクトル b と行列、テーブル、またはデータセット配列 x を受け入れます。この関数ハンドルは、x と同じ行数のベクトル f を返します。たとえば、関数ファイル hougen.m は次の計算を実行します。

hougen(b,x)=b(1)x(2)x(3)/b(5)1+b(2)x(1)+b(3)x(2)+b(4)x(3).

MATLAB® コマンド ラインに type hougen と入力して、関数を調べます。

function yhat = hougen(beta,x)
%HOUGEN Hougen-Watson model for reaction kinetics.
%   YHAT = HOUGEN(BETA,X) gives the predicted values of the
%   reaction rate, YHAT, as a function of the vector of 
%   parameters, BETA, and the matrix of data, X.
%   BETA must have 5 elements and X must have three
%   columns.
%
%   The model form is:
%   y = (b1*x2 - x3/b5)./(1+b2*x1+b3*x2+b4*x3)
%
%   Reference:
%      [1]  Bates, Douglas, and Watts, Donald, "Nonlinear
%      Regression Analysis and Its Applications", Wiley
%      1988 p. 271-272.

%   Copyright 1993-2004 The MathWorks, Inc. 
%   B.A. Jones 1-06-95.

b1 = beta(1);
b2 = beta(2);
b3 = beta(3);
b4 = beta(4);
b5 = beta(5);

x1 = x(:,1);
x2 = x(:,2);
x3 = x(:,3);

yhat = (b1*x2 - x3/b5)./(1+b2*x1+b3*x2+b4*x3);

hougen.m と同じ計算を実行する無名関数を作成することができます。

modelfun = @(b,x)(b(1)*x(:,2) - x(:,3)/b(5))./...
(1 + b(2)*x(:,1) + b(3)*x(:,2) + b(4)*x(:,3));

式の文字列表現

行列 X のデータとベクトル y の応答

  • X の最初の予測子 (列) として 'x1' を使用し、2 番目の予測子として 'x2' を使用するというように、式を表現します。

  • 最適化するパラメーターのベクトルを、'b1''b2' のように表現します。

  • 'y ~ (mathematical expressions)' のように、式を作成します。

たとえば、反応データへの応答は次のように表現します。

modelfun = 'y ~ (b1*x2 - x3/b5)/(1 + b2*x1 + b3*x2 + b4*x3)';

テーブルまたはデータセット配列のデータの場合、そのテーブルまたはデータセット配列の変数名をもつ文字列で表現した式を使用することができます。式の左側に応答変数名を配置し、その後に ~ を挿入し、さらに応答式を表す文字列を続けます。

次の例は、データセット配列内の reaction データに対する応答を表す文字列の作成方法を示しています。

  1. reaction データを読み込みます。

    load reaction
  2. データセット配列にデータを挿入します。ここで各変数に xn または yn の文字列で指定されている名前を付けます。

    ds = dataset({reactants,xn(1,:),xn(2,:),xn(3,:)},...
        {rate,yn});
  3. データセット配列の 1 行目を調べます。

    ds(1,:)
    
    ans = 
    
        Hydrogen    n_Pentane    Isopentane    ReactionRate
        470         300          10            8.55 
  4. データセット配列内の名前を使用して、hougen 式を作成します。

    modelfun = ['ReactionRate ~ (b1*n_Pentane - Isopentane/b5) /'...
    ' (1 + Hydrogen*b2 + n_Pentane*b3 + Isopentane*b4)']
    
    modelfun =
    ReactionRate ~ (b1*n_Pentane - Isopentane/b5) / ...
         (1 + Hydrogen*b2 + n_Pentane*b3 + Isopentane*b4)

初期ベクトル beta0 の選択

近似反復の初期ベクトルである beta0 は、結果の近似モデルの品質に大きな影響を与えることがあります。beta0 は問題の次元数を指定するため、正しい長さにする必要があります。beta0 を適切に選択すると、高速で信頼性の高いモデルが作成されますが、選択が不適切な場合は、計算に時間がかかるか、不適切なモデルが作成される可能性があります。

適切な beta0 の選択についてアドバイスするのは難しいでしょう。このベクトルの特定の成分が正または負であると考えられる場合は、beta0 にこれらの特性を設定します。他の成分の近似値が既知の場合は、それらの値を beta0 に含めます。ただし、適切な値がわからない場合は、次のような乱数ベクトルを試してください。

beta0 = randn(nVars,1);
% or
beta0 = 10*rand(nVars,1);

非線形モデルでデータを近似

テーブルまたはデータセット配列 tbl を使用して非線形回帰モデルで近似するには、次の構文を使用します。

mdl = fitnlm(tbl,modelfun,beta0)

データセット配列 X と数値応答ベクトル y を使用して非線形回帰モデルを近似するには、次の構文を使用します。

mdl = fitnlm(X,y,modelfun,beta0)

入力パラメーターの表現については、データの準備非線形モデルの表現および初期ベクトル beta0 の選択を参照してください。

fitnlm はテーブルまたはデータセット配列 tbl 内の応答変数が最後の列であると仮定しています。これを変更するには、ResponseVar 名前と値のペアを使用して応答列に名前を付けます。

品質の調査と近似非線形モデルの調整

いくつかの診断プロットは、モデルの品質を調べるのに役立ちます。plotDiagnostics(mdl) は、てこ比のプロットやクックの距離プロットなど、さまざまなプロットを提供します。plotResiduals(mdl) は近似モデルとそのデータの差異を示します。

mdl のプロパティも、モデルの品質に関係しています。mdl.RMSE はデータと近似モデルの間の二乗平均平方根誤差を示します。mdl.Residuals.Raw は生の残差を示します。mdl.Diagnostics に含まれるフィールドは、Leverage や CooksDistance など、特に興味深い観測値を特定するのに役立ちます。

次の例は、診断、残差およびスライス プロットを使用して近似線形モデルを調べる方法を示しています。

標本データを読み込みます。

load reaction

関数 hougen.m を使用して、rate の非線形モデルを reactants の関数として作成します。

beta0 = ones(5,1);
mdl = fitnlm(reactants,...
    rate,@hougen,beta0);

データとモデルのてこ比のプロットを作成します。

plotDiagnostics(mdl)

てこ比の高い点が 1 つあります。その点を特定します。

[~,maxl] = max(mdl.Diagnostics.Leverage)
maxl =

     6

残差プロットを調べます。

plotResiduals(mdl,'fitted')

顕著な外れ値はありません。

モデルの各予測子の効果を知るために、スライス プロットを使用します。

plotSlice(mdl)

青い縦の破線をドラッグすると、応答の 1 つの予測子を変更した効果を調べることができます。たとえば、X2 の線を右へドラッグすると、X3 の線の傾斜が変化することがわかります。

非線形モデルを使用した応答の予測またはシミュレーション

この例では、predictfevalrandom の各メソッドを使用して、新しいデータへの応答を予測してシミュレーションする方法を説明します。

コーシー分布からランダムに標本を生成します。

rng('default')
X = rand(100,1);
X = tan(pi*X - pi/2);

モデル y = b1*(pi /2 + atan((x � b2) / b3)) に従って応答を生成し、その応答にノイズを追加します。

modelfun = @(b,x) b(1) * ...
    (pi/2 + atan((x - b(2))/b(3)));
y = modelfun([12 5 10],X) + randn(100,1);

任意パラメーターの b = [1,1,1] から始まるモデルを近似します。

beta0 = [1 1 1]; % An arbitrary guess
mdl = fitnlm(X,y,modelfun,beta0)

mdl = 


Nonlinear regression model:
    y ~ b1*(pi/2 + atan((x - b2)/b3))

Estimated Coefficients:
          Estimate      SE       tStat       pValue  
          ________    _______    ______    __________

    b1    12.082      0.80028    15.097    3.3151e-27
    b2    5.0603       1.0825    4.6747    9.5063e-06
    b3      9.64      0.46499    20.732    2.0382e-37


Number of observations: 100, Error degrees of freedom: 97
Root Mean Squared Error: 1.02
R-Squared: 0.92,  Adjusted R-Squared 0.918
F-statistic vs. zero model: 6.45e+03, p-value = 1.72e-111

近似値はパラメーター [12,5,10] の数パーセント内に収まります。

近似を検証します。

plotSlice(mdl)

predict

メソッド predict は平均応答の予測と、要求された場合には信頼限界を提供します。X 値 [-15;5;12] での応答について予測された応答値と予測された信頼区間を求めます。

Xnew = [-15;5;12];
[ynew,ynewci] = predict(mdl,Xnew)
ynew =

    5.4122
   18.9022
   26.5161


ynewci =

    4.8233    6.0010
   18.4555   19.3490
   25.0170   28.0151

信頼区間はスライス プロットに反映されます。

feval

メソッド feval は平均応答を予測します。データセット配列からモデルを構成する場合、predict を使用するより feval を使用したほうが便利なことがよくあります。

データセット配列から非線形モデルを作成します。

ds = dataset({X,'X'},{y,'y'});
mdl2 = fitnlm(ds,modelfun,beta0);

X 値 [-15;5;12] における予測モデル応答 (CDF) を求めます。

Xnew = [-15;5;12];
ynew = feval(mdl2,Xnew)
ynew =

    5.4122
   18.9022
   26.5161

random

random メソッドは新しいランダム応答値をシミュレートします。これは、平均予測に学習データと同じ分散値のランダム外乱を加えたものと同じです。

Xnew = [-15;5;12];
ysim = random(mdl,Xnew)
ysim =

    6.0505
   19.0893
   25.4647

メソッド random を再実行します。結果が変わります。

ysim = random(mdl,Xnew)
ysim =

    6.3813
   19.2157
   26.6541

非線形回帰ワークフロー

この例では、一般的な非線形回帰ワークフローの使用方法を示します。データのインポート、非線形回帰、品質テスト、品質改善のための修正、そしてモデルに基づく予測という流れです。

手順 1. データを準備します。

データ reaction を読み込みます。

load reaction

ワークスペース内のデータを調べます。reactants は 13 行 3 列の行列です。各行は 1 つの観測に対応し、各列は 1 つの変数に対応します。変数名は次の xn で指定されています。

xn
xn =

Hydrogen  
n-Pentane 
Isopentane

同様に、rate は 13 の応答を含むベクトルで、変数名は次の yn で指定されています。

yn
yn =

Reaction Rate

ファイル hougen.m に含まれる反応速度の非線形モデルは、3 つの予測子変数をもつ関数です。5 次元ベクトル b と 3 次元ベクトル x の場合は、次のようになります。

hougen(b,x)=b(1)x(2)x(3)/b(5)1+b(2)x(1)+b(3)x(2)+b(4)x(3).

解の開始点として、いずれかのベクトルに b を使用します。

beta0 = ones(5,1);

手順 2. 非線形モデルをデータに近似する。

mdl = fitnlm(reactants,...
    rate,@hougen,beta0)
mdl = 

Nonlinear regression model:
    y ~ hougen(b,X)

Estimated Coefficients:
          Estimate    SE          tStat     pValue 
    b1      1.2526     0.86702    1.4447    0.18654
    b2    0.062776    0.043562    1.4411    0.18753
    b3    0.040048    0.030885    1.2967    0.23089
    b4     0.11242    0.075158    1.4957    0.17309
    b5      1.1914     0.83671    1.4239     0.1923

Number of observations: 13, Error degrees of freedom: 8
Root Mean Squared Error: 0.193
R-Squared: 0.999,  Adjusted R-Squared 0.998
F-statistic vs. constant model: 1.81e+03, p-value = 7.36e-12

手順 3. モデルの品質を調べる。

二乗平均平方根誤差は、観測値の範囲に比べて比較的小さい値です。

[mdl.RMSE min(rate) max(rate)]
ans =

    0.1933   0.0200   14.3900

残差プロットを調べます。

plotResiduals(mdl)

モデルはこのデータに適切なようです。

診断プロットを調べて、外れ値を探します。

plotDiagnostics(mdl,'cookd')

観測 6 はプロット ラインから外れているようです。

手順 4. 外れ値を削除する。

Exclude 名前と値のペアを使用して、近似から外れ値を除外します。

mdl1 = fitnlm(reactants,...
    rate,@hougen,ones(5,1),'Exclude',6)
mdl1 = 

Nonlinear regression model:
    y ~ hougen(b,X)

Estimated Coefficients:
          Estimate    SE          tStat     pValue 
    b1       0.619      0.4552    1.3598    0.21605
    b2    0.030377    0.023061    1.3172    0.22924
    b3    0.018927     0.01574    1.2024    0.26828
    b4    0.053411    0.041084       1.3    0.23476
    b5      2.4125      1.7903    1.3475     0.2198

Number of observations: 12, Error degrees of freedom: 7
Root Mean Squared Error: 0.198
R-Squared: 0.999,  Adjusted R-Squared 0.998
F-statistic vs. constant model: 1.24e+03, p-value = 4.73e-10

モデル係数が mdl の係数から大きく変わりました。

手順 5. 両方のモデルのスライス プロットを調べる。

応答に対する各予測子の効果を調べるには、plotSlice(mdl) を使用してスライス プロットを作成します。

plotSlice(mdl)

plotSlice(mdl1)

これらのプロットは非常によく似ていますが、mdl1 では信頼限界の幅がわずかに広くなります。近似しているデータ点が 1 つ少なく、観測が 7% 以上少ないことが示されるため、このような違いが生じることは理解できます。

手順 6. 新しいデータで予測する。

新しいデータを作成し、両方のモデルの応答を予測します。

Xnew =  [200,200,200;100,200,100;500,50,5];
[ypred yci] = predict(mdl,Xnew)
ypred =

    1.8762
    6.2793
    1.6718


yci =

    1.6283    2.1242
    5.9789    6.5797
    1.5589    1.7846
[ypred1 yci1] = predict(mdl1,Xnew)
ypred1 =

    1.8984
    6.2555
    1.6594


yci1 =

    1.6260    2.1708
    5.9323    6.5787
    1.5345    1.7843

モデル係数の類似度は低いにもかかわらず、予測はほぼ同じ結果となります。

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