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

Теоретическая часть

Метод максимального или наибольшего правдоподобия предложен Р. Фишером [6, 13]. С помощью этого метода производится точечная оценка неизвестных параметров априорно известного закона распределения случайной величины.

Рассмотрим сначала суть метода при оценке параметров дискретного распределения случайной величины [6].

Обозначим вероятность того, что в результате испытания величина $$X$$ примет значение $$x_{i}\mbox{ }(i=1,2,...,n)$$, через $$p_i(x_i,\theta)$$.

Определение. Функцией правдоподобия случайной дискретной величины $$X$$ называют функцию аргумента $$\theta$$:

где $$x_{1},\mbox{ }x_{2},...,x_{n}$$ — фиксированные числа, полученные при измерении случайной величины $$X$$.

В качестве точечной оценки параметра $$\theta$$ принимают такое его значение $$\theta^{*}=\theta^{*}(x_1,x_2,...,x_n)$$, при котором функция правдоподобия достигает максимума. Оценку $$\theta^{*}$$ называют оценкой максимального правдоподобия.

Для упрощения расчетов в рассмотрение вводится логарифм функции правдоподобия $$lnL$$, которую называют логарифмической функцией правдоподобия. Функции $$L$$ и $$lnL$$ достигают максимума при одном и том же значении своего аргумента, поэтому вместо отыскания максимума функции $$L$$ ищут максимум функции $$lnL$$. Записывая необходимое условие экстремума функции правдоподобия в случае скалярного параметра, получаем уравнения правдоподобия

или

где $$\vec x_n$$ — заданная выборка случайных величин.

Уравнение правдоподобия (7.3) с логарифмической функцией, как правило, более простое относительно функции правдоподобия (7.2).

Если распределение случайной величины $$Х$$ зависит от вектора параметров $$\vec \theta=(\theta_1,\theta_2,...,\theta_r)$$, то уравнение (7.3) заменяется системой уравнений

Именно уравнения (7.3) и (7.4) принято называть уравнениями правдоподобия [13]. Во многих случаях решение системы (7.4), являющейся, как правило, нелинейной, приходится искать численными методами.

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

Пусть $$X$$ — непрерывная случайная величина, которая в результате $$n$$ испытаний приняла значения $$x_{1},\mbox{ }x_{2},...,x_{n}$$. Предполагается, что вид плотности распределения $$f(x)$$ задан, но неизвестен параметр $$\theta$$, которым определяется эта функция.

Определение. Функцией правдоподобия непрерывной случайной величины $$Х$$ называют функцию аргумента $$\theta$$

где $$x_{1},\mbox{ }x_{2},...,x_{n}$$ — фиксированные числа.

Оценку максимального правдоподобия неизвестного параметра $$\theta$$ распределения непрерывной случайной величины ищут так же, как в случае дискретной величины.

Замечание. Если плотность распределения непрерывной случайной величины $$Х$$ определяется двумя неизвестными параметрами $$\theta_1$$ и $$\theta_2$$, то функция правдоподобия является функцией двух независимых аргументов $$\theta_1$$ и $$\theta_2$$:

Как для дискретных распределений, так и для непрерывных точку максимума логарифмической функции распределения $$lnL$$ аргумента $$\theta$$ можно искать через необходимое условие экстремума:

Найденную точку максимума $$\theta^{*}$$ принимают в качестве оценки максимального правдоподобия параметра $$\theta$$.

Метод максимального правдоподобия имеет ряд достоинств: его оценки, вообще говоря, состоятельны (но они могут быть смещенными), распределены асимптотически нормально (при больших значениях $$n$$ приближенно нормально) и имеют наименьшую дисперсию по сравнению с другими асимптотически нормальными оценками; если для оцениваемого параметра $$\theta$$ существует эффективная оценка $$\theta^{*}$$, то уравнение правдоподобия имеет единственное решение $$\theta^{*}$$ ; этот метод наиболее полно использует данные выборки об оцениваемом параметре, поэтому он особенно полезен в случае малых выборок. Недостаток метода состоит в том, что он часто требует сложных вычислений.

Практическая часть

1. Оценка параметра экспоненциального распределения

Рассматривается пример поиска методом максимального правдоподобия оценки параметра $$\lambda$$ экспоненциального распределения случайной величины, для которой функция плотности имеет вид

К характеристикам экспоненциального распределения относятся математическое ожидание $$M[X]$$ и дисперсия $$D[X]$$:

Замечание. Во встроенных функциях MATLAB параметром экспоненциального распределения является математическое ожидание случайной величины.

Возможная программная реализация точечной оценки параметра экспоненциального распределения:

clear,clc,close all
%%% Проверка на закрытие диалоговых окон
try
   global h11
   close(h11);
end
try
   global n11
   
close(n11);
end
 
try
   global v11
   close(v11)
end
 
%% ВВОД ТЕОРЕТИЧЕСКОГО ПАРАМЕТРА РАСПРЕДЕЛЕНИЯ
options.Resize = 'on';
options.WindowStyle = 'modal'; %%'normal';
options.Interpreter = 'tex';
P1 = inputdlg({'\bfВвод параметра:......................................................'},...
sprintf('Теоретическая величина параметра'),1,{'1.23'},options);
%% ПРЕОБРАЗОВАНИЕ К СТРОКОВОЙ ПЕРЕМЕННОЙ
P2 = char(P1);
%% ПРЕОБРАЗОВАНИЕ К ЧИСЛУ С ДВОЙНОЙ ТОЧНОСТЬЮ
P0 = str2num(P2);
 
%% КОНТРОЛЬ ВВОДА ПАРАМЕТРА
if isempty(P0)
h11 = errordlg('Параметр должен быть действительным положительным числом!','Ошибка ввода');
    return
end
%% КОНТРОЛЬ ВВОДА ПАРАМЕТРА
global h11
if P0 <= 0 | ~isreal(P0) | ~isfinite(P0)
    h11 = errordlg('Параметр должен быть конечным действительным положительным числом!','Ошибка ввода');
    return
end
% ВВОД ЧИСЛА ПРОГОНОВ ПРОГРАММЫ
n1 = inputdlg({'\bfВвод числа прогонов программы..........................'},...
    'Число прогонов программы',1,{'10'}, options);
 
% ПРЕОБРАЗОВАНИЕ К ЧИСЛОВОЙ ПЕРЕМЕННОЙ
n = str2num(char(n1));

%% Контроль ввода цифр
if isempty(n)
    global n11
n11 = errordlg('Число прогонов программы должно быть целым положительным числом!', 'Ошибка ввода');    
    return
end
if ~isreal(n) | ~isfinite(n)
global n11

n11 = errordlg('Число прогонов программы должно быть целым положительным числом!', 'Ошибка ввода');   
    return
end
%% Контроль целого положительного числа циклов
if n <= 0 | n ~= round(n)
global n11
n11 = errordlg('Число прогонов программы должно быть целым положительным числом!', 'Ошибка ввода');
return
end
 
% ВВОД ЧИСЛА ИЗМЕРЕНИЙ СЛУЧАЙНОЙ ВЕЛИЧИНЫ
v1 = inputdlg({'\bfВвод числа измерений случайной величины...................................'},...
    'Число измерений случайной величины',1,{'1234'}, options);
 
% ПРЕОБРАЗОВАНИЕ К ЧИСЛОВОЙ ПЕРЕМЕННОЙ
v = str2num(char(v1));
if isempty(v)
    global v11
v11 = errordlg('Число измерений должно быть положительным целым числом!','Ошибка ввода');   
    return
end
if ~isreal(v) | ~isfinite(v)
  global v11
v11 = errordlg('Число измерений должно быть положительным целым числом!','Ошибка ввода');   
    return
end
% КОНТРОЛЬ ЦЕЛОГО ЧИСЛА ИЗМЕРЕНИЙ 
% СЛУЧАЙНОЙ ВЕЛИЧИНЫ
if v <= 0 | v ~= round(v)
    global v11
v11 = errordlg('Число измерений должно быть положительным целым числом!','Ошибка ввода');
return
        end
syms m
k = 0;
%% ЦИКЛ ЗАДАННОГО ЧИСЛА ПРОГОНОВ ПРОГРАММЫ
for I = 1:n
    k=k+1;
%% ФОРМИРОВАНИЕ ЧИСЛА ИЗМЕРЕНИЙ СЛУЧАЙНОЙ ВЕЛИЧИНЫ
t = exprnd(1/P0,v,1);
%% ФОРМИРОВАНИЕ ФУНКЦИИ МАКСИМАЛЬНОГО 
%% ПРАВДОПОДОБИЯ
L = m^(length(t))*exp(-m*sum(t));
%% ЛОГАРИФМИЧЕСКАЯ ФУНКЦИЯ МАКСИМАЛЬНОГО 
%% ПРАВДОПОДОБИЯ
Lg = log(L);

%% ДИФФЕРЕНЦИРОВАНИЕ
dLg = diff(Lg,m);

%% ПРЕОБРАЗОВАНИЕ СИМВОЛЬНОЙ ПЕРЕМЕННОЙ К СТРОКОВОЙ 
dLg = char(dLg);
%% РЕШЕНИЕ УРАВНЕНИЯ ОТНОСИТЕЛЬНО ОЦЕНИВАЕМОГО 
%% ПАРАМЕТРА
as1(k) = double(solve(dLg));
%% УСРЕДНЕНИЕ ОЦЕНИВАЕМОГО ПАРАМЕТРА
as(k) = mean(as1);
end
%% ОКОНЧАНИЕ ЦИКЛА ЗАДАННОГО ЧИСЛА ПРОГОНОВ ПРОГРАММЫ
mcp = mean(as);
%% ВЫВОД РЕЗУЛЬТАТОВ В КОМАНДНОЕ ОКНО
fprintf('\n\t%s%g\n \t%s%g\n','Теоретический параметр: ',P0,...
'Оценка параметра: ', mcp)
fprintf('\tОтносительная погрешность: %g%s\n',abs(P0-mcp)/P0*100,'%')
 
%% ГРАФИЧЕСКИЕ ПОСТРОЕНИЯ
figure(1)
%% set(gcf,'position',[a(3)/90,a(3)/20,a(3)/2.1,a(4)/2])
plot(1:n,as1,'r:','linew',2),grid off,hold on,
plot(1:n,as,'linew',2),
title(sprintf('%s%g','\bfТеоретический параметр\fontsize{12} \lambda\fontsize{10} = ',P0))
 xlabel('\bf Количество циклов'),
ylabel('\bf Эмпирический параметр\fontsize{14} \lambda'),
legend('\bf Измеряемая величина\fontsize{12} \lambda',...
    '\bf Средняя величина\fontsize{12} \lambda'),
set(gcf,'color','w')
 
%% ПОСТРОЕНИЕ ТЕОРЕТИЧЕСКОЙ И ЭМПИРИЧЕСКОЙ 
%% ФУНКЦИИ ПЛОТНОСТИ
t = 0 : 0.1 : 4;
y1 = P0*exp(-P0*t); %exppdf(t,1/P0); % встроенная функция
y2 = mcp*exp(-mcp*t); %exppdf(t,1/mcp);
 
figure(2)
plot(t, y1, 'r', 'linew',2),
hold on
plot(t, y2, 'bo', 'linew',2)
grid off
legend('\bf Теоретическая функция плотности (PDF)',...
'\bf Эмпирическая функция плотности'),
text(t(end)/3,2/3*max(max([y1,y2])),['\bf',...
sprintf('Теоретический параметр: %g\n Эмпирический параметр: %g',P0,mcp)])
xlabel('\bf Случайная величина'),
ylabel('\bf Функция плотности'),
set(gcf,'color','w')

Задание 1

№ 1: n = 200; № 2: n = 300; № 3: n = 400; № 4: n = 500; № 5: n = 600;
№ 6: n = 700; № 7: n = 800; № 8: n = 900; № 9: n = 1000; № 10: n = 2000;

№ 1: (10-19); № 2: (20-29); № 3: (30-39); № 4: (40-49); № 5 (50-59);
№ 6: (60-69); № 7: (70-79); № 8: (80-89); № 9: (90-99); № 10: (100-110).