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

Открыть 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'

The 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. The axes with title Histogram of residuals contains an object of type patch.

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

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

plotDiagnostics(mdl,'cookd')

Figure contains an axes. The axes 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 and other objects of type uimenu, uicontrol. Axes 1 contains 5 objects of type line. Axes 2 contains 5 objects of type line. Axes 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

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

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