Данные и модель для подгонки

Мы будем использовать данные, собранные для изучения загрязнения воды промышленными и бытовыми отходами. Эти данные подробно описаны в Box, G.P., W.G. Hunter и J.S. Hunter, Statistics for Experimenters (Wiley, 1978, pp. 483-487). Ответная переменная представляет собой биохимическую потребность в кислороде в мг/л, а предикторная переменная представляет собой время инкубации в днях.

x = [1 2 3 5 7 10]';
y = [109 149 149 191 213 224]';

plot(x,y,'ko');
xlabel('Incubation (days), x'); ylabel('Biochemical oxygen demand (mg/l), y');

Мы предполагаем, что известно, что первые два наблюдения были сделаны с меньшей точностью, чем остальные наблюдения. Например, они могли быть изготовлены с помощью другого инструмента. Другой распространенной причиной взвешивания данных является то, что каждое зарегистрированное наблюдение является фактически средним из нескольких измерений, выполненных при одном и том же значении x. В данных здесь предположим, что первые два значения представляют одно необработанное измерение, в то время как остальные четыре представляют собой среднее из 5 необработанных измерений. Тогда было бы уместно взвешивать по количеству измерений, которые шли в каждое наблюдение.

w = [1 1 5 5 5 5]';

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

modelFun = @(b,x) b(1).*(1-exp(-b(2).*x));

Только основываясь на грубой визуальной посадке, кажется, что кривая, проведенная через точки, может выровняться на уровне около 240 где-то в районе x = 15. Поэтому мы будем использовать 240 в качестве начального значения для b1, и поскольку e ^ (-.5 * 15) мало по сравнению с 1, мы будем использовать .5 в качестве начального значения для b2.

start = [240; .5];

Подгонка модели без весов

Опасность игнорирования ошибки измерения заключается в том, что на посадку могут оказывать чрезмерное влияние неточные измерения и, следовательно, могут не обеспечить хорошую подгонку к измерениям, которые точно известны. Давайте поместим данные без весов и сравним их с точками.

nlm = fitnlm(x,y,modelFun,start);
xx = linspace(0,12)';
line(xx,predict(nlm,xx),'linestyle','--','color','k')

Подгонка модели с весами

Обратите внимание, что подогнанная кривая тянется к первым двум точкам, но, похоже, упускает тренд других точек. Попробуем повторить посадку с помощью весов.

wnlm = fitnlm(x,y,modelFun,start,'Weight',w)
line(xx,predict(wnlm,xx),'color','b')

wnlm = 


Nonlinear regression model:
    y ~ b1*(1 - exp( - b2*x))

Estimated Coefficients:
          Estimate       SE       tStat       pValue  
          ________    ________    ______    __________

    b1     225.17         10.7    21.045    3.0134e-05
    b2    0.40078     0.064296    6.2333     0.0033745


Number of observations: 6, Error degrees of freedom: 4
Root Mean Squared Error: 24
R-Squared: 0.908,  Adjusted R-Squared 0.885
F-statistic vs. zero model: 696, p-value = 8.2e-06

Оценочное стандартное отклонение совокупности в этом случае описывает среднее изменение для «стандартного» наблюдения с весом или точностью измерения 1.

wnlm.RMSE

ans =

   24.0096

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

coefCI(wnlm)

ans =

  195.4650  254.8788
    0.2223    0.5793

Оценка кривой отклика

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

[ypred,ypredci] = predict(wnlm,xx,'Simultaneous',true);
plot(x,y,'ko', xx,ypred,'b-', xx,ypredci,'b:');
xlabel('x'); ylabel('y');
legend({'Data', 'Weighted fit', '95% Confidence Limits'},'location','SouthEast');

Обратите внимание, что две точки с пониженным весом не совпадают по кривой с остальными точками. Это так, как вы ожидали для взвешенной посадки.

Также можно оценить интервалы прогнозирования для будущих наблюдений при заданных значениях x. Эти интервалы будут фактически принимать вес или точность измерения, равную 1.

[ypred,ypredci] = predict(wnlm,xx,'Simultaneous',true,'Prediction','observation');
plot(x,y,'ko', xx,ypred,'b-', xx,ypredci,'b:');
xlabel('x'); ylabel('y');
legend({'Data', 'Weighted fit', '95% Prediction Limits'},'location','SouthEast');

Абсолютная шкала весов фактически не влияет на оценки параметров. Масштабирование весов на любую константу дало бы нам те же оценки. Но они влияют на доверительные границы, поскольку границы представляют собой наблюдение с весом 1. Здесь видно, что точки с большим весом кажутся слишком близкими к подогнанной линии, по сравнению с доверительными пределами.

Пока predict метод не позволяет нам изменять веса, мы можем сделать некоторую постобработку и исследовать, как кривая будет искать более точную оценку. Предположим, мы заинтересованы в новом наблюдении, которое основано на среднем из пяти измерений, как и последние четыре точки на этом графике. Мы могли бы уменьшить ширину интервалов в коэффициент sqrt (5).

halfwidth = ypredci(:,2)-ypred;
newwidth = halfwidth/sqrt(5);
newci = [ypred-newwidth, ypred+newwidth];
plot(x,y,'ko', xx,ypred,'b-', xx,newci,'r:');
xlabel('x'); ylabel('y');
legend({'Data', 'Weighted fit', 'Limits for weight=5'},'location','SouthEast');

Остаточный анализ

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

r = wnlm.Residuals.Raw;
plot(x,r.*sqrt(w),'b^');
xlabel('x'); ylabel('Residuals, yFit - y');

Есть некоторые свидетельства систематических закономерностей в этом остаточном сюжете. Обратите внимание, как последние четыре остатка имеют линейный тренд, предполагающий, что модель может увеличиваться недостаточно быстро по мере увеличения x. Кроме того, величина остатков имеет тенденцию уменьшаться по мере увеличения x, предполагая, что ошибка измерения может зависеть от x. Они заслуживают исследования, однако, есть так мало точек данных, что трудно придавать значение этим кажущимся шаблонам.