Выбор между фитингом кривой и распределительным фитингом

Фитинг кривой и фитинг распределения - это различные типы анализа данных.

Фитинг кривой

В следующих экспериментальных данных переменная предиктора time, время после приема препарата. Переменная ответа: concконцентрация препарата в кровотоке. Предположим, что только данные ответа conc влияет экспериментальная ошибка.

time = [ 0.1   0.1   0.3   0.3   1.3   1.7   2.1   2.6   3.9   3.9 ...
         5.1   5.6   6.2   6.4   7.7   8.1   8.2   8.9   9.0   9.5 ...
         9.6  10.2  10.3  10.8  11.2  11.2  11.2  11.7  12.1  12.3 ...
        12.3  13.1  13.2  13.4  13.7  14.0  14.3  15.4  16.1  16.1 ...
        16.4  16.4  16.7  16.7  17.5  17.6  18.1  18.5  19.3  19.7]';
conc = [0.01  0.08  0.13  0.16  0.55  0.90  1.11  1.62  1.79  1.59 ...
        1.83  1.68  2.09  2.17  2.66  2.08  2.26  1.65  1.70  2.39 ...
        2.08  2.02  1.65  1.96  1.91  1.30  1.62  1.57  1.32  1.56 ...
        1.36  1.05  1.29  1.32  1.20  1.10  0.88  0.63  0.69  0.69 ...
        0.49  0.53  0.42  0.48  0.41  0.27  0.36  0.33  0.17  0.20]';

Предположим, вы хотите смоделировать концентрацию крови как функцию времени. График conc против time.

plot(time,conc,'o');
xlabel('Time');
ylabel('Blood Concentration');

Предположим, что conc следует двухпараметрической кривой Вейбулла как функция time. Кривая Вейбулла имеет форму и параметры

где $$a$$ - горизонтальное масштабирование, $$b$$ - параметр формы, а $$c$$ - вертикальное масштабирование.

Подгонка модели Вейбулла с использованием нелинейных наименьших квадратов.

modelFun =  @(p,x) p(3) .* (x./p(1)).^(p(2)-1) .* exp(-(x./p(1)).^p(2));
startingVals = [10 2 5];
nlModel = fitnlm(time,conc,modelFun,startingVals);

Постройте график кривой Вейбулла на данные.

xgrid = linspace(0,20,100)';
line(xgrid,predict(nlModel,xgrid),'Color','r');

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

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

nlModel2 = fitnlm(time,log(conc),@(p,x) log(modelFun(p,x)),startingVals);

Добавьте новую кривую к существующему графику.

line(xgrid,exp(predict(nlModel2,xgrid)),'Color',[0 .5 0],'LineStyle','--');
legend({'Raw Data','Additive Errors Model','Multiplicative Errors Model'});

Объект модели nlModel2 содержит оценки точности. Лучшая практика - проверить благость модели. Например, создайте остаточные графики в логарифмической шкале, чтобы проверить допущение постоянной дисперсии на наличие мультипликативных ошибок.

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

Функции для фитинга кривой

Распределительный фитинг

Предположим, что требуется смоделировать распределение срока службы электрических компонентов. Переменная life измеряет время до отказа 50 идентичных электрических компонентов.

life = [ 6.2 16.1 16.3 19.0 12.2  8.1  8.8  5.9  7.3  8.2 ...
        16.1 12.8  9.8 11.3  5.1 10.8  6.7  1.2  8.3  2.3 ...
         4.3  2.9 14.8  4.6  3.1 13.6 14.5  5.2  5.7  6.5 ...
         5.3  6.4  3.5 11.4  9.3 12.4 18.3 15.9  4.0 10.4 ...
         8.7  3.0 12.1  3.9  6.5  3.4  8.5  0.9  9.9  7.9]';

Визуализация данных с помощью гистограммы.

binWidth = 2;
lastVal = ceil(max(life));
binEdges = 0:binWidth:lastVal+1;
h = histogram(life,binEdges);
xlabel('Time to Failure');
ylabel('Frequency');
ylim([0 10]);

Поскольку данные о времени жизни часто следуют распределению Вейбулла, одним из подходов может быть использование кривой Вейбулла из предыдущего примера подгонки кривой для подгонки гистограммы. Чтобы попробовать этот подход, преобразуйте гистограмму в набор точек (x, y), где x - центр ячейки, а y - высота ячейки, а затем поместите кривую в эти точки.

counts = histcounts(life,binEdges);
binCtrs = binEdges(1:end-1) + binWidth/2;
h.FaceColor = [.9 .9 .9];
hold on
plot(binCtrs,counts,'o');
hold off

Однако подгонка кривой к гистограмме является проблематичной и обычно не рекомендуется.

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

Однако $$\mathrm{b/a}$$ заменяет параметр масштаба $$c$$, поскольку функция должна интегрироваться в 1. Чтобы подогнать распределение Вейбулла к данным с использованием максимального правдоподобия, используйте fitdist и указать 'Weibull' в качестве имени дистрибутива. В отличие от наименьших квадратов, максимальное правдоподобие находит pdf-файл Вейбулла, который лучше всего соответствует масштабированной гистограмме без минимизации суммы квадратных разностей между высотой pdf-файла и высоты полосы.

pd = fitdist(life,'Weibull');

Постройте график масштабированной гистограммы данных и наложите соответствующий pdf.

h = histogram(life,binEdges,'Normalization','pdf','FaceColor',[.9 .9 .9]);
xlabel('Time to Failure');
ylabel('Probability Density');
ylim([0 0.1]);
xgrid = linspace(0,20,100)';
pdfEst = pdf(pd,xgrid);
line(xgrid,pdfEst)

Лучшая практика - проверить благость модели.

Хотя подгонка кривой к гистограмме обычно не рекомендуется, в некоторых случаях процесс подходит. Пример см. в разделе Подгонка пользовательских одномерных распределений.

Функции распределительного фитинга