非線形回帰

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

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

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 には表示されません。

非線形モデルの表現

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

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

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

無名関数または関数ファイルの関数ハンドル
式のテキスト表現

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

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

$hougen (b, x) = \frac{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 データに対する応答を表す文字ベクトルの作成方法を示しています。

reaction データを読み込みます。
```
load reaction
```
データセット配列にデータを格納します。各変数の名前は xn または yn で指定されているものになります。
```
ds = dataset({reactants,xn(1,:),xn(2,:),xn(3,:)},...
    {rate,yn});
```

データセット配列の 1 行目を調べます。

ds(1,:)

ans = 

    Hydrogen    n_Pentane    Isopentane    ReactionRate
    470         300          10            8.55

データセット配列内の名前を使用して、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)

Figure contains an axes object. The axes object with title Case order plot of leverage, xlabel Row number, ylabel Leverage contains 2 objects of type line. One or more of the lines displays its values using only markers These objects represent Leverage, Reference Line.

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

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

maxl = 
6

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

plotResiduals(mdl,'fitted')

Figure contains an axes object. The axes object with title Plot of residuals vs. fitted values, xlabel Fitted values, ylabel Residuals contains 2 objects of type line. One or more of the lines displays its values using only markers

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

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

plotSlice(mdl)

Figure Prediction Slice Plots contains 3 axes objects and another object of type uigridlayout. Axes object 1 contains 4 objects of type line, patch. Axes object 2 contains 4 objects of type line, patch. Axes object 3 contains 4 objects of type line, patch.

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

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

ライブスクリプトを開く

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

コーシー分布から無作為に標本を生成します。

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)

Figure Prediction Slice Plots contains an axes object and another object of type uigridlayout. The axes object contains 5 objects of type line, patch. One or more of the lines displays its values using only markers

predict

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

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

ynew = 3×1

    5.4122
   18.9022
   26.5161

ynewci = 3×2

    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 = 3×1

    5.4122
   18.9022
   26.5161

random

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

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

ysim = 3×1

    6.0505
   19.0893
   25.4647

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

ysim = random(mdl,Xnew)

ysim = 3×1

    6.3813
   19.2157
   26.6541