Задание 1. Построение оценки максимального правдоподобия и её доверительного интервала для параметра экспоненциального распределения длительности события

Экспоненциальное распределение характерно для длительности электронных приборов и элементов при отсутствии дрейфа параметров – старения. Длительность события, имеющая экспоненциальное распределение, описывается плотностью вероятности вида

и функцией дожития

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

, (3.1)

где , если наблюдение нецензурировано, , если наблюдение цензурировано, m – число нецензурированных наблюдений (случаев). Необходимым условием максимума выражения (3.1) является равенство нулю первой производной LnL по a

,

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

.

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

(3.2)

Убедимся, что найден именно максимум логарифма правдоподобия. Вычислим вторую производную LnL по a

.

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

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

1. При отсутствии цензурирования, построение доверительного интервала для оценки максимального правдоподобия параметра a основывается на асимптотической нормальности первой производной LnL по a. Согласно этому свойству, при большом числе наблюдений случайная величина при истинном значении параметра a имеет приближённо нормальное распределение с математическим ожиданием 0 и дисперсией . Вторая производная также вычисляется при истинном значении параметра a. Нормальность распределения обозначается записью . По свойству нормального распределения, при большом числе наблюдений для случайной величины , вычисленной при истинном значении параметра a, с вероятностью 0.95 справедливо

Выражение для первой производной логарифма правдоподобия без цензурирования запишем в виде

,

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

.

Подставляя найденные выражения в (3.3), получим

.

Преобразуя, получим, что с вероятностью 0.95 истинное значение параметра a лежит между границами

. (3.3)

2. При наличии цензурирования для построения доверительного интервала оценки параметра экспоненциального распределения можно использовать метод профиля логарифма правдоподобия. Этот метод основан на том, что при большом числе наблюдений случайная величина имеет распределение с одной степенью свободы. Из свойств распределения следует, что с вероятностью 0.95 справедливо неравенство .

Метод профиля логарифма правдоподобия иллюстрирует рис. 1. Вычисления проводятся в следующем порядке:

1. строится график функции ;

2. максимальное значение графика достигается при ;

3. вычисляется критическое значение и решается неравенство ;

4. интервал, соответствующий решению неравенства, с доверительной вероятностью 0.95 является доверительным интервалом для параметра экспоненциального распределения.

Рисунок 1. Иллюстрация вычисления доверительного интервала по профилю логарифма правдоподобия.

1.1 Получите у преподавателя значения параметров a и n. Запишите оператор для генерации независимой выборки значений экспоненциально распределённой случайной величины с параметром a

X=exprnd(1/a,1,n);

1.2 По формуле (3.2) вычислить оценку максимального правдоподобия для параметра a, обозначив её a0. Вывести на экран точное значение параметра а, полученное в задании, и построенную оценку a0. Для вывода использовать оператор

disp([‘a=’,num2str(a),’ estimate a=’,num2str(a0)]);

1.3 По формуле (3.3) вычислить границы 95% доверительного интервала для параметра a, обозначив их a1 и a2. Вывести a1 и a2 на экран компьютера с помощью оператора

disp([num2str(a1),’ < a < ‘,num2str(a2)]);

1.4 Считая, что цензурирования нет, построить по формуле (3.1) график логарифма правдоподобия, вычислить критическое значение t, найти приближённое решение неравенства , вывести результат на экран компьютера. Для вычислений можно использовать следующие операторы

at1=0.05; at2=0.2; % величины at1 и at2 могут изменяться

% в зависимости от сгенерированной выборки

at=at1:0.001:at2;

i=0;

for a=at

i=i+1;

LnL(i)=n*log(a)-a*sum(X);

end

t=n*log(a0)-a0*sum(X)-1.92;

plot(at,LnL,'k',[at(1),at(end)],[t t],’--’)

box off; xlabel('a')

at=at(LnL>=t);

disp([num2str(at(1)),’ < a < ‘,num2str(at(end))]);

1.5 С помощью команды run Lab3 выполните записанные операторы. Значения параметров, построенный график и выведенные на экран числа внесите в отчёт о выполнении лабораторной работы. Сделайте вывод о близости оценки параметра экспоненциального распределения, построенной методом максимального правдоподобия, к истинному значению параметра. Убедитесь, что истинное значение параметра экспоненциального распределения попадает в 95% доверительный интервал, вычисленный двумя методами.

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

Модель Вейбулла используется при описании длительности безотказной работы технических систем при наличии постепенных отказов, вызванных дрейфом параметров системы в процессе эксплуатации. Случайная величина, имеющая распределение, которое задаётся моделью Вэйбулла, описывается функцией распределения

,

где a – масштабный множитель, k – параметр формы. Плотность вероятности имеет вид

.

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

, (3.4)

где , если наблюдение нецензурировано, , если наблюдение цензурировано, m – число нецензурированных наблюдений (случаев). Необходимым условием максимума выражения (3.4) является равенство нулю первых частных производных LnL по a и по k

,

откуда получаем, что оценка максимального правдоподобия параметра масштаба a равна

(3.5).

Для получения оценки максимального правдоподобия параметра формы k необходимо решить другое уравнение. Запишем

.

Подставив выражение для оценки , получим уравнение для определения оценки k, которое решается численными методами

(3.6).

Подставляя найденное значение оценки в (3.5) получим оценку параметра .

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

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

Метод профиля логарифма правдоподобия иллюстрирует рис. 3.2. Вычисления проводятся в следующем порядке:

1. строится поверхность функции ;

2. максимальное значение поверхности достигается при ;

3. вычисляется критическое значение и определяется область, задаваемая неравенством ;

4. найденная область является 95% доверительной областью для параметров модели Вейбулла.

Рисунок 3.2. Иллюстрация вычисления доверительной области по поверхности логарифма правдоподобия.

2.1 Получите у преподавателя значения параметров a, k и n. Запишите оператор для генерации независимой выборки значений случайной величины, имеющей распределение, описываемое моделью Вейбулла с параметрами a и k

X=wblrnd(1/a,k,1,n);

2.2 Считая, что цензурирования нет, построить по формуле (3.4) поверхность логарифма правдоподобия,

% значения переменных a1, a2, k1 и k2 могут зависеть от сгенерированных данных и параметров модели Вейбулла

a1=1.5;a2=3; k1=0.8; k2=1.4;

at=a1:0.01:a2;

kt=k1:0.01:k2;

LnL=zeros(length(at),length(kt));

for i=1:length(at)

for j=1:length(kt)

LnL(i,j)=n*log(at(i)*kt(j))+(kt(j)-1)*sum(log(at(i)*X))-sum((at(i)*X).^k(j));

end

end

2.3 Вычислить критическое значение t, нарисовать область, ограничивающую множество решений неравенства

maxLnL=max(LnL(:));

t=maxLnL-3.0;

contour(a,k,LnL',[lev lev]);

xlabel(‘a’); ylabel(‘k’);

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

ah=a(max(LnL,[],2)==maxLnL); % оценка масштаба

kh=k(max(LnL,[],1)==maxLnL); % оценка параметра формы

disp([‘a=’,num2str(ah),’ k=’,num2str(kh)])

2.5 Нанести на плоскости, построенной в пункте 2.3, символ * зелёного цвета с координатами a и k, а также символ О красного цвета с координатами ah и kh

hold on;

plot(a,k,’*g’);

plot(ah,kh,’Or’);

2.6 С помощью команды run Lab3 выполните записанные операторы. Значения параметров, построенный график и выведенные на экран числа внесите в отчёт о выполнении лабораторной работы. Сделайте вывод о близости оценок параметров модели Вейбулла, построенных методом максимального правдоподобия, к их истинным значениям.