Нелинейный рабочий процесс регрессии

Скрипт Open Live Script

В этом примере показано, как сделать типичный нелинейный рабочий процесс регрессии: импортируйте данные, соответствуйте нелинейной регрессии, протестируйте ее качество, измените их, чтобы улучшить качество и сделать предсказания на основе модели.

Шаг 1. Подготовьте данные.

Загрузите reaction данные.

load reaction

Исследуйте данные в рабочей области. reactants матрица с 13 строками и 3 столбцами. Каждая строка соответствует одному наблюдению, и каждый столбец соответствует одной переменной. Имена переменных находятся в xn:

xn

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

Точно так же rate вектор из 13 ответов, с именем переменной в yn:

yn

yn = 
'Reaction Rate'

hougen.m файл содержит нелинейную модель скорости реакции в зависимости от этих трех переменных предикторов. Для 5-D вектора $b$ и 3-D вектор $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 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 contains 2 objects of type line. 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. Это различие понятно, поскольку существует тот меньше точки данных в подгонке, представляя на более чем 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

Даже при том, что коэффициенты модели отличаются, предсказания почти идентичны.

Документация