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

Мы будем использовать данные, собранные для изучения загрязнения воды, вызванного промышленными и бытовыми отходами. Эти данные подробно описаны в 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. Они заслуживают исследования, однако существует так мало точек данных, что трудно придать значимость этим очевидным шаблонам.