"Компьютерные технологии в медико-биологических системах"
Тема: «Оценка параметров модели Гомпертца-Мейкхема методом максимального правдоподобия»
Оглавление
Введение2
Глава 1. Основные положения математической теории описания функции дожития3
Глава 2. Методы описания дожития8
2.1. Метод максимального правдоподобия как классический метод математической статистики8
2.2. Описание модели Гомпертца-Мейкхема9
2.3. Оценка параметров модели Гомпертца-Мейкхема методом максимального правдоподобия10
2.4. Метод перепараметризации функции11
Глава 3. Построение модели Гомпертца-Мейкхема по данным дожития с помощью программы MATLAB12
Заключение21
Список литературы22
Введение
Начиная со второй половины ХХ века, наблюдается постоянное увеличение продолжительности жизни населения экономически развитых стран. Среди стран с наибольшей средней продолжительностью жизни, средняя продолжительность жизни за четыре года увеличивается на один год. Для выяснения причин такого роста и построения надёжного прогноза ожидаемой в будущем продолжительности жизни ученые изучают статистику смертности населения. Задача ученых – с помощью имеющейся статистики, найти закономерности и законы, по которым возможно приблизительно рассчитать вероятность дожития или смерти отдельно взятого человека, находящегося при заданных условиях проживания.
Существует множество различных моделей описания дожития и смертности. Каждая из них справедлива при определенных условиях, на каком-то определенном возрастном промежутке и описывается набором параметров. Среди этих моделей популярны модели экспоненциального дожития, модель Вейбулла и модель Гомпертца-Мейкхама.
В курсовой работе рассматривается оценивание параметров в модели Гомпертца-Мейкхема. Рассматриваются методы приближения распределения с неизвестными параметрами к распределению случайной выборки. С помощью программы MATLAB реализована программа-функция, максимально приближающая функции дожития и распределения модели Гомпертца-Мейкхема к распределению имеющейся случайной выборки методом максимального правдоподобия.
Приводится теоретическое обоснование применяемых методов и графическое представление результатов.
Глава 1. Основные положения математической теории описания функции дожития
Функция дожития и распределения
Формальное рассмотрение дожития[1].
Пусть T – длительность процесса (случайной величины)
- функция дожития.
- функция распределения
Тогда из этих понятий выходят и их определения.
Функция дожития это вероятность того, что какой-либо объект доживет до момента времени .
Функция распределения это вероятность того, что какой-либо объект не доживет до момента времени .
Равномерное распределение
Экспоненциальное распределение
Плотность распределения вероятности наступления события
Плотность распределения вероятности – вероятность попасть в интервал.
Интенсивность наступления события
Интенсивность наступления события или мгновенный риск определяется уравнениями:
Кумулятивный риск
Из второго определения мгновенного риска можно выразить функцию дожития .
В полученном выражении подставив пределы интегрирования можно ввести еще одно понятие.
Интеграл, стоящий после минуса в экспоненте называется кумулятивным риском :
Тогда выходит еще одно определение понятия функции дожития:
Модели дожития
|
Название
|
|
|
|
|
|
Экспоненциальное распределение |
||||
|
Модель Гомпертца |
||||
|
Модель Вейбулла |
||||
|
Логарифмически-нормальное распределение |
|
|||
|
Модель Гомпертца-Мейкхема |
Где
– начальная смертность
– темп старения
– внешний риск
параметр масштаба
параметр формы
Цензурирование
Цензурирование процесса – это прерывание наблюдения процесса, в связи с неизвестным временем его окончания.
Цензурирование – прекращение наблюдения до наступления события.
Пусть
T – длительность процесса
C – длительность интервала наблюдения
T, C – независимы.
Глава 2. Методы описания дожития
2.1. Метод максимального правдоподобия как классический метод математической статистики
Предположим, что группой ученых был проведен один и тот же опыт или эксперимент в одинаковых условиях n раз.
В каждом подобном опыте было зафиксировано свое значение случайной величины. Все значения независимы друг от друга.
Предполагается, что данная выборка подчиняется априори известному закону с неизвестным параметром .
Задача – оценить значение параметра .
Допустим, что функция плотности распределения данного процесса есть
, где - параметр
Тогда, для подбора нужного нам значения параметра , вводится функция, которая называется функцией правдоподобия и которая определяется выражением
Смысл ее таков, что при нахождении максимума данной функции, мы определим параметр , при котором априорный закон будет совпадать с законом имеющейся выборки.
– оценка максимального правдоподобия.
Обычно, для нахождения максимального значения функции правдоподобия, берут ее натуральный логарифм.
При переходе от самой функции правдоподобия к ее логарифму максимальное значение, в простых случаях остается в той же точке. Поэтому данный переход вполне справедлив и при этом, в большинстве случаев, упрощается нахождение максимального значения параметра .
В простых случаях максимум вычисляется аналитически.
Например, имеется выборка , подчиняющаяся экспоненциальному распределению. Тогда
Выразив через мгновенный риск и функцию дожития, функция правдоподобия примет следующий вид:
Возьмем натуральный логарифм от полученного выражения:
Подставив мгновенный и кумулятивный риски экспоненциального распределения в получившееся уравнение, получим:
Максимум функции правдоподобия находится при помощи первой и второй производных.
Одной производной не достаточно в связи с тем, что функция может иметь не один экстремум. Поэтому если вторая производная в точке экстремума меньше нуля, то она является локальным максимумом.
Приравниваем первую производную функции правдоподобия к нулю:
Отсюда
- оценка максимального правдоподобия
Проверим полученную точку экстремума с помощью второй производной:
Следовательно - значение параметра, при котором функция правдоподобия имеет максимальное значение.
2.2. Описание модели Гомпертца-Мейкхема
Закон смертности Гомпертца — Мейкхама (иногда просто Закон Гомпертца, Распределение Гомпертца) — статистическое распределение, которое описывает смертность человека и большинства многоплодных животных. Согласно закону Гомпертца — Мейкхама, смертность является суммой независимого от возраста компонента (члена Мейкхама) и компонента, зависимого от возраста (функция Гомпертца), который экспоненциально возрастает с возрастом и описывает старение организма. В защищённых средах, где внешние причины смерти отсутствуют (в лабораторных условиях, в зоопарках или для людей в развитых странах) независимый от возраста компонент часто становится малым, и формула упрощается до функции Гомпертца. Распределение было получено и опубликовано актуарием и математиком Бенджамином Гомпертцем в 1832 году[2].
Согласно закону Гомпертца — Мейкхама, мгновенный риск смерти за фиксованный короткий промежуток времени после достижения возраста x составляет:
где t — возраст, a, b, c — коэффициенты. Таким образом, размер популяции снижается с возрастом по удвоенной экспоненте:
Закон смертности Гомпертца — Мейкхама наилучшим образом описывает динамику смертности человека в диапазоне возраста 30—80 лет. В области большего возраста смертность не возрастает так быстро, как предусматривается этим законом смертности.
Исторически смертность человека до 1950-х годов была в большей мере вызвана независимым от времени компонентом закона смертности (членом или параметром Мейкхама), тогда как зависимый от возраста компонент (функция Гомпертца) почти не изменялась. После 1950-х годов картина изменилась, что привело к снижению смертности в позднем возрасте и так называемой «де-ректангуляризации» (сглаживанию) кривой выживания.
2.3. Оценка параметров модели Гомпертца-Мейкхема методом максимального правдоподобия
Имеется случайная выборка xi, имеющая распределение согласно модели Гомпертца-Мейкхама. Требуется определить значения параметров a, b, c, при которых распределение Гомпертца-Мейкхама будет максимально подобным распределению случайной выборки[1].
Имеем цензурированную выборку xi, имеющую распределение согласно модели Гомпертца-Мейкхема:
|
цензурировано |
Мгновенный риск и функция дожития модели Гомпертца-Мейкхема определяются по формулам:
Где a, b, c – параметры, x – возраст.
Значения параметров a, b, c определяются методом максимального правдоподобия, при котором требуется достичь максимизации функции правдоподобия. При введении этого понятия в формуле также учитываются цензурированные выборки.
Где – плотность распределения вероятности, зависящей от параметра и выборки xi
Как известно, плотность вероятности определяется по формуле:
Тогда функция правдоподобия будет иметь вид:
При взятии натурального логарифма от функции правдоподобия положение максимума не изменится. Тогда:
Где – кумулятивный риск.
2.4. Метод перепараметризации функции
При численном нахождении параметров с помощью функции правдоподобия, могут получаться отрицательные значения.
В теории дожития все параметры функций лежат в положительной области.
Для устранения данной проблемы используется метод перепараметризации, при котором параметр заменяют экспонентой какого-то нового параметра. Например, для распределения Гомпертца-Мейкхема:
Тогда максимизация функции правдоподобия проводится по параметрам p1, p2, p3, которые выражаются из исходных параметров a, b, c.
Тогда при любом значении параметров p1, p2, p3 значения исходных параметров всегда будут положительными.
Глава 3. Построение модели Гомпертца-Мейкхема по данным дожития с помощью программы MATLAB
Задание
Создать программу-функцию на MATLAB для вычисления и графического представления функции дожития и смертности для модели Гомпертца-Мейкхема с параметрами, рассчитанными по данным дожития.
Вариант 1
- a=0.1, b=0, c=0.001, n=25
- a=0.1, b=0, c=0.5, n=150
- a=0.1, b=0, c=0.001, n=1000
Напишем программу-функцию, которая будет:
а) генерировать случайную выборку
б) цензурировать случайную выборку
в) максимизировать функцию правдоподобия
г) выводить на экран графики функций дожития и смертности, заданных исходными параметрами и их оценками
Так как мы имеем дело с программой-функцией, ввод параметров осуществляется через её формальные параметры.
То есть, чтобы вызвать программу-функцию, например, с первыми значениями исходных параметров, в командном окне требуется ввести:
kurs(0.1, 0, 0.001, 25, 5)
Код программы:
% 1) Ввод исходных параметров a0, b0, c0, n и m выполняются через формальные параметры программы-функции
function kurs(a0, b0, c0, n, m)
%a0 - начальная смертность
%b0 - темп старения
%с0 - внешний риск
%n - количество элементов в выборке
%m - количество цензурированных элементов выборки
% 2) Генерирование случайной независимой выборки объёма n длительностей, имеющих распределение согласно модели Гомпертца-Мейкхема
q=1;
x=exprnd(1/q,n,1);
%q - среднее значение выборки
%x - результирующая выборка
% 3) Цензурирование
aux=randperm(n);
z=aux(1:m);
В следующем пункте программы речь пойдет о максимизации натурального логарифма от функции правдоподобия.
В MATLAB'е эта операция выполняется с помощью нахождения минимума минуса функции командой fminsearch.
Так как при взятии функции со знаком «минус» исходная функция зеркально отразится относительно оси х, модуль значения её в точке локального экстремума не изменится и данная операция будет вполне уместной для достижения наших целей. А именно, определение значения не самого экстремума, а точки, в которой он находится.
% 4) Минимализация функциии -lnL
[p2, F]=fminsearch(@f, log([a0, b0, c0]), [], x, z);
% 5) Вывод результатов
a=exp(p2(1))
b=exp(p2(2))
c=exp(p2(3))
% 6) Графический вывод функций
% укажем промежуток, на котором удобно рассмотреть графики
x=0:50;
% 6.1 Графики функций дожития
В ходе выполнения программы параметр b может оказаться нулевым или очень малым (<=1.0e-5). Так как параметр находится в знаменателе, подобные значения нельзя подставлять в функцию дожития S1(t) и распределения F1(t).
Та же проблема может возникнуть, если мы изначально зададим очень маленький или нулевой параметр b0. Тогда функции дожития S0(t) и распределения F0(t) также не будут иметь смысл.
Чтобы устранить данную проблему создаются условия, что, если параметр b и b0 меньше какого-то предельного значения, функции дожития и распределения будут вычисляться по экспоненциальной модели.
% функция дожития S0(t) для исходных значений параметров
if b0<=1.0e-5
S0=exp(-a0*x-c0*x);
else
S0=exp(-a0*(exp(b0*x)-1)/b0-c0*x);
end
% функция дожития S1(t) для оценочных значений параметров
if b<=1.0e-5
S1=exp(-a*x-c*x);
else
S1=exp(-a*(exp(b*x)-1)/b-c*x);
end
plot(x, S0, 'r-');
hold on;
plot(x, S1, 'k--');
title('Функции дожития','Fontname','Arial cyr');
xlabel('t');
ylabel('S(t)');
legend('S0(t)','S1(t)');
% 6.2 Графики функций распределения
figure;
% функция распределения для исходных значений параметров
if b0<=1.0e-5
F0=1-exp(-a0*x-c0*x);
else
F0=1-exp(-a0*(exp(b0*x)-1)/b0-c0*x);
end
% функция распределения для оценочных значений параметров
if b<=1.0e-5
F1=1-exp(-a*x-c*x);
else
F1=1-exp(-a*(exp(b*x)-1)/b-c*x);
end
plot(x, F0, 'r-');
hold on;
plot(x, F1, 'k--');
title('Функции смертности','Fontname','Arial cyr');
xlabel('t');
ylabel('F(t)');
legend('F0(t)','F1(t)');
Код программы-функции f:
%Функция правдоподобия для модели Гомпертца-Мейкхема
function y=f(p, x, z)
% p - параметр
% x - случайная выборка, полученная в программе-функции Kurs.m
% z - цензурированная выборка, полученная в программе-функции Kurs.m
Можно заметить, что здесь применен метод перепараметризации. В функции правдоподобия параметры a, b, c заменены на экспоненты параметров p(1), p(2), p(3). Тогда минимизация функции должна вестись по натуральному логарифму исходных параметров, что и делается в 4-ом пункте программы-функции Kurs.m.
a=exp(p(1)); b=exp(p(2)); c=exp(p(3));
x1=x;
x1(z)=[];
d=exp(b*x1);
r=exp(b*x);
y=sum(log(a*d)+c);
if b<1.0e-5
y=y-sum(a*x+c*x);
else
y=y-sum(a*(r-1)/b+c*x);
end
y=-y;
Графический вывод результатов:
Рис.1 Графики функций дожития первого набора значений
(S0(t) – функция дожития исходных параметров, S1(t) – функция дожития оценок параметров)
Рис.2 Графики функций распределения первого набора значений
(F0(t) – функция дожития исходных параметров, F1(t) – функция дожития оценок параметров)
Рис.3 Графики функций дожития второго набора значений
(S0(t) – функция дожития исходных параметров, S1(t) – функция дожития оценок параметров)
Рис.4 Графики функций распределения второго набора значений
(F0(t) – функция дожития исходных параметров, F1(t) – функция дожития оценок параметров)
Рис.5 Графики функций дожития третьего набора значений
(S0(t) – функция дожития исходных параметров, S1(t) – функция дожития оценок параметров)
Рис.6 Графики функций распределения третьего набора значений
(F0(t) – функция дожития исходных параметров, F1(t) – функция дожития оценок параметров)
В достаточно близком приближении графиков замечается их неточность в совпадение. Полное совпадение графиков возможно только тогда, когда функция правдоподобия стремиться к бесконечности.
Заключение
В ходе выполнения курсовой работы были проведены:
- Полный обзор основных положений математической теории описания функции дожития
- Изучение основных методов описания дожития
- Реализация программы-функции в MATLAB, максимально приближающую функцию дожития и распределения модели Гомпертца-Мейкхема к распределению имеющейся случайной выборки методом максимального правдоподобия
- Графический вывод, построенных программой MATLAB, функций дожития и распределения модели Гомпертца-Мейкхема
При каждом выводе результатов программы-функции Kurs.m значения параметров модели Гомпертца-Мейкхема были близки со значениями исходных параметров.
Также графический вывод показывает, что функции дожития и распределения оценок параметров близки к тем же функциям случайной выборки.
Все это доказывает правильность выполнения практической части курсовой работы.
Список литературы:
- Лекции по курсу: «Компьютерные технологии в медико-биологических системах», 2014 год
- Интернет-ресурс:
http://ru.wikipedia.org/wiki/Закон_смертности_Гомпертца-Макегама, 2014год
Скачать курсовую работу: