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

Мы будем использовать данные, собранные, чтобы изучить загрязнение воды, вызванное промышленными и бытовыми отходами. Эти данные описаны подробно в Поле, G.P., В.Г. Хантер, и Дж.С. Хантер, Статистика для Экспериментаторов (Вайли, 1978, стр 483-487). Переменная отклика является биохимическим спросом на кислород в mg/l, и переменный предиктор является временем инкубации в днях.

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

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

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

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