1.4.3. Обоснование закона распределения

Осуществим выбор закона распределения между следующими типами подходящих по теоретически соображениям распределений, обычно [43] применяемыми для решения подобных задач:

• логарифмическое,

• логнормальное,

• гамма-распределение,

• распределение Вейбулла,

• экспоненциальное распределение,

• распределение Рэллея,

• распределение Гомпертца,

• логлогистическое.

Функция распределения (с точностью до параметров) может быть известна из теоретических соображений. В таком случае задача вычисления параметров распределений может трактоваться как идентификация математической модели. Параметрическая статистическая модель выживания описывается с помощью следующих функций:

• функция плотности вероятности распределения

• кумулятивная или интегральная функция распределения - функция распределения длительностей до момента отказа t,

• функция выживания – вероятность безотказной работы до момента t,

• функция интенсивности отказов (функция риска) .

Под длительностью t может пониматься время работоспособности («жизни»), количество циклов до отказа и т.п., в зависимости от конкретной задачи.

Теоретические значения частот (непрерывных) распределений элементарно вычисляются как произведение плотности теоретического распределения на численность выборки и на длину соответствующего классового интервала.

Общая методика

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

Реализация метода начинается с составления функции максимального правдоподобия (ФМП) в виде

где - вектор неизвестных параметров статистической модели,

- массив длительностей,

, - соответствующий массив индикаторов цензурирования (для конкретной длительности индикаторов равен 1, если выборка цензурирована),

n – численность массива длительностей.

Отметим, что в дальнейших выкладках число цензурированных длительностей стандартно обозначено как

Вектор искомых параметров находится из условия максимума ФМП:

Максимум ФМП находится из условия равенства нулю частных производных ФМП по искомым параметрам, т.е. искомые параметры удовлетворяют уравнениям

Для упрощения максимизируют не саму ФМП, а логарифм ФМП. Эта возможность основана на том факте, что ФМП и логарифм достигают максимума при одних и тех же значениях искомых параметров, однако работать с логарифмом ФМП значительно проще:

Задача сводится, таким образом, к аналитическому либо численному (одним из методов оптимизации) решению полученной линейной или нелинейной системы уравнений.

Логарифмические модели

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

,

где массив логарифмов длительностей,

– вектор параметров,

параметр положения,

параметр масштаба.

Логарифмическая ФМП может быть записана как

где массив стандартизованных длительностей.

В дальнейших выкладках понадобятся следующие очевидные выражения для производных

Тогда компоненты вектора градиента логарифмической ФМП параметрам

Компоненты матрицы вторых производных логарифмической ФМП по параметрам (матрицы Гессе) вычисляются как

С учетом введенных обозначений интерационная схема максимизации логарифмической ФМП алгоритма метода Ньютона-Рафсона может быть записана как

номер интерации.

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

Логнормальное распределение

Плотность логнормального распределения (логарифмически нормального) распределения [78] с двумя параметрами имеет вид

Введем нормированную величину

Тогда можно записать

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

Соответствующая функция выживания

где Ф(.) – функция стандартного нормального распределения.

С учетом нормировки функция максимального правдоподобия (ФМП) запишется как

Соответствующая логарифмическая ФМП имеет вид

Производные, необходимые для интерационной схемы метода Ньютона-Рафсона, запишутся как

Для моделирования обычно используется связь с нормальным распределением. Поэтому, достаточно сгенерировать нормально распределённую случайную величину, но в случае анализа выживаемости не всегда удается сгенерировать нормально распределённую случайную величину и вычислить ее экспоненту.

Гамма-распределение

Плотность гамма-распределение имеет вид

Соответствующая функция выживания

где I(.,.) – неполная гамма-функция.

Поэтому функция максимального правдоподобия (ФМП) запишется как

Соответствующая логарифмичеcкая ФМП имеет вид

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

Распределение Вейбулла

Плотность распределение Вейбулла с двумя параметрами имеет вид

Соответствующая функция выживания

Поэтому функция максимального правдоподобия (ФМП) запишется как

После преобразований окончательно получаем

Соответствующая логарифмическая ФМП имеет вид

Логарифмическая ФМП примет окончательный вид

Для вычисления значений искомых параметров найдем частные производные логарифмической ФМП по искомым параметрам и приравняем к нулю. Сначала найдем производную по :

Отсюда уравнение для вычисления параметра получается как

Вычислив производную по параметру ,

с учетом выражения для параметра , получаем нелинейное уравнение для поиска параметра методом деления отрезка пополам.

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

Экпоненциальное распределение

Плотность экспоненциального распределения имеет вид

Поэтому функция максимального правдоподобия (ФМП) запишется как

После преобразований окончательно получаем

Соответствующая логарифмическая ФМП имеет вид

Для вычисления значения искомого параметра найдем производную логарифмической ФМП по данному параметру и приравняем ее к нулю. Производную по параметру имеет вид:

Отсюда уравнение для вычисления параметра получается как

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

Распределение Рэллея

Плотность распределения Рэлея имеет вид

Соответствующая функция выживания

Поэтому функция максимального правдоподобия запишется как

После преобразований окончательно получаем

Соответствующая логарифмическая ФМП имеет вид

Для вычисления значения искомого параметра найдем производную логарифмической ФМП по данному параметру и приравняем ее к нулю. Производная по параметру имеет вид:

Уравнение для вычисления параметра получается как

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

Распределение Гомпертца

Плотность распределения Гомпертца имеет вид

Соответствующая функция выживания

Поэтому максимальная функция правдоподобия запишется как

После преобразований получаем

Соответствующая логарифмическая ФМП имеет вид

lnL

Для вычислений значений искомых параметров найдем частные производные логарифмической ФМП по искомым параметрам и приравняем их нулю. Сначала найдем производную по :

Отсюда уравнение для вычисления параметра получается как

Вычислив производную по параметру , получаем

С учетом выражения для параметра , получаем нелинейное уравнение для поиска параметра в виде:

Решение уравнения может быть произведено одним из методов оптимизации – в простейшем случае методом деления отрезка пополам.

Логлогистическое распределение

Плотность логлогистического распределения (логарифмически логистического) распределения с двумя параметрами имеет вид

Введем нормированную величину

Тогда можно записать

Соответствующая функция выживания

С учетом нормировки функция максимального правдоподобия (ФМП) запишется как

Соответствующая логарифмическая ФМП имеет вид

Производные, необходимые для итерационной метода Ньютона-Рафсона, запишутся как

Как было сказано ранее, значения частот (непрерывных) распределений элементарно вычисляются как произведение плотности теоретического распределения на численность выборки и на длину соответствующего классового интервала. Сделаем подбор распределения длительностей (рис 2.1).

Выбор распределения длительностей Логнормальное распределение ( ) Гамма-распределение ( , K=1,2401) Распределение Вейбулла 0,0166 1,213193) Экспоненциальное распределение (гамма) 0,032765 Распределение Рэлея (бета) 20,35534 Распределение Гомпертца (альфа, бета) -0,10566 0,0882 Логлогистическое распределение (мю,сигма) 3,032074 0,666967

Таблица 1– Выбор распределения длительности

Класс	Опыт	Логнормал	Гамма	Вейбулл	Экспоненц	Рэлей	Гомпертц	Логлогистич
5	17	21,23787	16,53556	16,1254	17,80108	7,493611	23,6279	18,17609
13	24	15,82921	14,84986	15,31512	13,69644	16,37559	7,662278	16,36322
21	14	10,54192	11,89591	12,64805	10,53826	19,05124	2,916585	11,42243
29	14	7,274097	9,177464	9,84838	8,108309	16,23566	1,189248	7,770953
37	7	5,219198	6,946953	7,394086	6,238665	10,95367	0,499461	5,401072
45	3	3,869969	5,198447	5,407352	4,80013	6,036344	0,212442	3,871878
53	1	2,947902	3,86015	3,87376	3,693298	2,760315	0,090855	2,861379
Хи-квт		14,50695	11,60679	10,09174	15,93746	21,33379	347,0922	11,09548
P-знач		0,005841	0,020528	0,03891	0,007025	0,0007	0	0,055512
AIC		66,12838	64,34876	63,86297	63,21139	61,00983	91,01506	61,07089

Рисунок 4 - Подбор распределения длительности

Достаточно важным вопросом является соответствие результатов расчета объективными критериями согласия для каждого распределения. Адекватной является модель с Р – значением, большим 0,05. В этом случае теоретическое и эмпирическое распределения значимо не различаются. Качество статистической модели можно оценить (также сравнить между собой различные модели), используя информационный критерий Акаике [101]

где - оценка логарифма функции максимального правдоподобия,

- вектор оценок параметров статистической модели,

K – число параметров модели.

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

где , - эмпирический массив длительности,

, – соответствующий массив индикаторов цензурирования,

n – численность массива длительностей,

- оценка функции плотности вероятности распределения,

- оценка соответствующей функции выживания.

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

Исходя из описаний всех видов распределений можно увидеть, что у каждого распределения есть недостатки, по причине которых их использование в риск-анализе будет не эффективно, либо не возможно (низкое быстродействие примитивных методов расчетов, сложность выполнения аналитического представления производных логарифмической ФМП). Адекватной является модель с Р – значением, большим 0,05. В этом случае теоретическое и эмпирическое распределения значимо не различаются. Качество статистической модели можно оценить (также сравнить между собой различные модели), используя информационный критерий Акаике. При сравнении нескольких статистических моделей лучшей считается модель с наименьшим значением AIC. Из рассчитанных параметров (таблица 2.1) видно, что самое подходящее распределение – логлогистическое, так как у него Р – значение = 0,055512, что больше 0,05 и AIC=61,07089, наименьший из представленных. В соответствие с построенным подбором распределений (рисунок 2.1), можно сделать предположение о том что, наиболее подходящим распределением является логлогистическое.