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

Открыть сценарий в реальном времени

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

Шаг 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$ ,

$hougen (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

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

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