非線形回帰ワークフロー

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

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

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

load reaction

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

xn

xn = 3x10 char array
    'Hydrogen  '
    'n-Pentane '
    'Isopentane'

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

yn

yn = 
'Reaction Rate'

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

$h o u g e n (b, x) = \frac{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. zero model: 3.91e+03, p-value = 2.54e-13

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

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

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

ans = 1×3

    0.1933    0.0200   14.3900

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

plotResiduals(mdl)

Figure contains an axes object. The axes object with title Histogram of residuals, xlabel Residuals, ylabel Probability density contains an object of type patch.

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

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

plotDiagnostics(mdl,'cookd')

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

観測 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. zero model: 2.67e+03, p-value = 2.54e-11

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

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

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

plotSlice(mdl)

Figure Prediction Slice Plots contains 3 axes objects and other objects of type uimenu, uicontrol. Axes object 1 contains 5 objects of type line. Axes object 2 contains 5 objects of type line. Axes object 3 contains 5 objects of type line.

plotSlice(mdl1)

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

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

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

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

ypred = 3×1

    1.8762
    6.2793
    1.6718

yci = 3×2

    1.6283    2.1242
    5.9789    6.5797
    1.5589    1.7846

[ypred1 yci1] = predict(mdl1,Xnew)

ypred1 = 3×1

    1.8984
    6.2555
    1.6594

yci1 = 3×2

    1.6260    2.1708
    5.9323    6.5787
    1.5345    1.7843

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