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

Скрипт 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)

Модель кажется достаточной для данных.

Исследуйте диагностический график искать выбросы.

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

Коэффициенты модели, измененные вполне немного от тех в mdl.

Шаг 5. Исследуйте графики среза обеих моделей.

Чтобы видеть эффект каждого предиктора на ответе, сделайте график среза с помощью plotSlice(mdl).

plotSlice(mdl)

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

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

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