非線形回帰
パラメトリック非線形回帰モデルとは
パラメトリック非線形モデルは、連続応答変数と複数の連続予測子変数との関係を次の形式で表します。
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 の応答データか、table tbl 内の入力データと tbl の列としての応答データで始まります。入力データの各行が、1 つの観測値を表します。各列が 1 つの予測子 (変数) を表します。
table tbl では、次のように 'ResponseVar' の名前と値のペアで応答変数を示します。
mdl = fitlm(tbl,'ResponseVar','BloodPressure');
応答変数は既定で最後の列です。
非線形回帰の場合、"カテゴリカル" 予測子は使用できません。カテゴリカル予測子は可能性のある固定セットからの値をとります。
入力データと応答データの両方で、欠損値を NaN で表します。
入力および応答データの table
Excel® スプレッドシートから table を作成するには、次のようにします。
tbl = readtable('hospital.xls',... 'ReadRowNames',true);
ワークスペース変数から table を作成するには、次のようにします。
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 は次の計算を実行します。
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)';
table 内のデータの場合、その table の変数名で表現した式を使用することができます。式の左側に応答変数名を配置し、その後に ~ を挿入し、さらに応答式を表す文字ベクトルを続けます。
次の例は、table 内の reaction データに対する応答を表す文字ベクトルの作成方法を示しています。
reactionデータを読み込みます。load reaction
table にデータを格納します。変数の名前は
xnおよびynにある名前に基づきます。tbl = table(reactants(:,1),reactants(:,2),reactants(:,3),... rate, VariableNames=["Hydrogen","n_Pentane",... "Isopentane","ReactionRate"])table を調べます。
tbl tbl = Hydrogen n_Pentane Isopentane ReactionRate ________ _________ __________ ____________ 470 300 10 8.55 285 80 10 3.79 470 300 120 4.82 470 80 120 0.02 470 80 10 2.75 100 190 10 14.39 100 80 65 2.54 470 190 65 4.35 100 300 54 13 100 300 120 8.5table 内の名前を使用して、
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);
非線形モデルのデータへの当てはめ
table tbl を使用して非線形回帰モデルを当てはめるには、次の構文を使用します。
mdl = fitnlm(tbl,modelfun,beta0)
データセット配列 X と数値応答ベクトル y を使用して非線形回帰モデルを当てはめるには、次の構文を使用します。
mdl = fitnlm(X,y,modelfun,beta0)
入力パラメーターの表現については、データの準備、非線形モデルの表現および初期ベクトル beta0 の選択を参照してください。
fitnlm は table 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 の線の傾斜が変化することがわかります。
非線形モデルを使用した応答の予測またはシミュレーション
この例では、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)
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 メソッドは平均応答の予測を提供します。
table から非線形モデルを作成します。
tbl = table(X,y,VariableNames=["X","Y"]); mdl2 = fitnlm(tbl,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
