Оценки наибольшего правдоподобия
phat = mle(data,Name,Value)
Например, можно задать тип распределения при помощи одних из этих аргументов name-value: Distribution, pdf, logpdf, или nloglf.
Чтобы вычислить MLEs для встроенного распределения, задайте тип распределения при помощи Distribution. Например, 'Distribution'\beta
Чтобы вычислить MLEs для пользовательского дистрибутива, задайте распределение при помощи pdf, logpdf, или nloglf, и задайте начальные значения параметров при помощи Start.
Найдите MLEs для встроенного распределения, что вы задаете использование Distribution аргумент значения имени.
Загрузите выборочные данные.
load carbigПеременная MPG содержит мили на галлон для различных моделей автомобилей.
Чертите гистограмму MPG данные.
histogram(MPG)
Распределение несколько правильно скошенный. Симметричное распределение, такое как нормальное распределение, не может быть подходящим вариантом.
Оцените параметры распределения Типа XII Шума для MPG данные.
phat = mle(MPG,'Distribution','burr')
phat = 1×3
34.6447 3.7898 3.5722
MLE для масштабного коэффициента α 34.6447. Оценки для двух параметров формы и из Шума распределение Типа XII 3.7898 и 3.5722, соответственно.
Сгенерируйте 100 случайных наблюдений от биномиального распределения с количеством испытаний = 20 и вероятность успеха = 0.75.
rng('default') % For reproducibility data = binornd(20,0.75,100,1);
Оцените вероятность успеха и 99%-х пределов достоверности с помощью симулированных выборочных данных. Необходимо задать количество испытаний (NTrials) для биномиального распределения.
[phat,pci] = mle(data,'Distribution','binomial','NTrials',20, ... 'Alpha',.01)
phat = 0.7615
pci = 2×1
0.7361
0.7856
Оценка вероятности успеха 0.7615, и нижние и верхние пределы 99%-го доверительного интервала 0.7361 и 0.7856, соответственно. Этот интервал покрывает истинное значение, используемое, чтобы симулировать данные.
Сгенерируйте выборочные данные размера 1000 от нецентрального распределения хи-квадрат со степенями свободы 8 и параметр нецентрированности 3.
rng default % for reproducibility x = ncx2rnd(8,3,1000,1);
Оцените параметры нецентрального распределения хи-квадрат от выборочных данных. Distribution аргумент значения имени не поддерживает нецентральное распределение хи-квадрат. Поэтому необходимо задать пользовательский нецентральный хи-квадрат PDF с помощью pdf аргумент значения имени и ncx2pdf функция. Необходимо также задать начальные значения параметров (Start аргумент значения имени) для пользовательского дистрибутива.
[phat,pci] = mle(x,'pdf',@(x,v,d)ncx2pdf(x,v,d),'Start',[1,1])
phat = 1×2
8.1052 2.6693
pci = 2×2
7.1120 1.6025
9.0983 3.7362
Оценка для степеней свободы 8.1052, и параметр нецентрированности 2.6693. 95%-й доверительный интервал для степеней свободы (7.1120,9.0983), и интервал для параметра нецентрированности (1.6025,3.7362). Доверительные интервалы включают истинные значения параметров 8 и 3, соответственно.
Загрузите выборочные данные.
load('readmissiontimes.mat');Данные включают ReadmissionTime, который имеет времена повторного доступа для 100 пациентов. Эти данные симулированы.
Задайте пользовательский журнал PDF для распределения Weibull с масштабным коэффициентом lambda и параметр формы k.
custlogpdf = @(data,lambda,k) ...
log(k) - k*log(lambda) + (k-1)*log(data) - (data/lambda).^k;Оцените параметры пользовательского дистрибутива и задайте его начальные значения параметров (Start аргумент значения имени).
phat = mle(ReadmissionTime,'logpdf',custlogpdf,'Start',[1,0.75])
phat = 1×2
7.5727 1.4540
Шкала и параметры формы пользовательского дистрибутива 7.5727 и 1.4540, соответственно.
Загрузите выборочные данные.
load('readmissiontimes.mat')Данные включают ReadmissionTime, который имеет времена повторного доступа для 100 пациентов. Эти данные симулированы.
Задайте пользовательскую отрицательную функцию логарифмической правдоподобности для распределения Пуассона параметром lambda, где 1/lambda среднее значение распределения. Необходимо задать функцию, чтобы принять логический вектор из информации о цензуре и целочисленный вектор из частот данных, даже если вы не используете эти значения в пользовательской функции.
custnloglf = @(lambda,data,cens,freq) ... - length(data)*log(lambda) + sum(lambda*data,'omitnan');
Оцените параметр пользовательского дистрибутива и задайте его начальное значение параметров (Start аргумент значения имени).
phat = mle(ReadmissionTime,'nloglf',custnloglf,'Start',0.05)
phat = 0.1462
Сгенерируйте выборочные данные размера 1000 от нецентрального распределения хи-квадрат со степенями свободы 10 и параметр нецентрированности 5.
rng('default') % For reproducibility x = ncx2rnd(10,5,1000,1);
Предположим, что параметр нецентрированности фиксируется в значении 5. Оцените степени свободы нецентрального распределения хи-квадрат от выборочных данных. Для этого задайте пользовательский нецентральный хи-квадрат PDF с помощью pdf аргумент значения имени.
[phat,pci] = mle(x,'pdf',@(x,v)ncx2pdf(x,v,5),'Start',1)
phat = 9.9307
pci = 2×1
9.5626
10.2989
Оценка для параметра нецентрированности 9.9307, и нижние и верхние пределы 95%-го доверительного интервала 9.5626 и 10.2989. Доверительный интервал включает истинное значение параметров 10.
Добавьте масштабный коэффициент в распределение хи-квадрат для адаптации к шкале данных и соответствуйте распределению.
Сгенерируйте выборочные данные размера 1000 от распределения хи-квадрат со степенями свободы 5 и масштабируйте данные на коэффициент 100.
rng default % For reproducibility x = 100*chi2rnd(5,1000,1);
Оцените степени свободы и масштабный коэффициент. Для этого задайте пользовательскую функцию плотности вероятности хи-квадрата с помощью pdf аргумент значения имени. Функция плотности требует a фактор для данных, масштабируемых .
[phat,pci] = mle(x,'pdf',@(x,v,s)chi2pdf(x/s,v)/s,'Start',[1,200])
phat = 1×2
5.1079 99.1681
pci = 2×2
4.6862 90.1215
5.5297 108.2146
Оценка для степеней свободы 5.1079, и шкала 99.1681. 95%-й доверительный интервал для степеней свободы (4.6862,5.5279), и интервал для масштабного коэффициента (90.1215,108.2146). Доверительные интервалы включают истинные значения параметров 5 и 100, соответственно.
Загрузите выборочные данные.
load('readmissiontimes.mat');Данные включают ReadmissionTime, который имеет времена повторного доступа для 100 пациентов. Вектор-столбец Censored содержит информацию о цензуре для каждого пациента, где 1 указывает, что подвергнутое цензуре правом наблюдение, и 0 указывает, что точное время повторного доступа наблюдается. Эти данные симулированы.
Задайте пользовательскую функцию плотности вероятности (PDF) и кумулятивная функция распределения (cdf) для экспоненциального распределения параметром lambda, где 1/lambda среднее значение распределения. Чтобы соответствовать распределению к подвергнутому цензуре набору данных, необходимо передать и PDF и cdf к mle функция.
custpdf = @(data,lambda) lambda*exp(-lambda*data); custcdf = @(data,lambda) 1-exp(-lambda*data);
Оцените параметр lambda из пользовательского дистрибутива для подвергнутых цензуре выборочных данных. Задайте начальное значение параметров (Start аргумент значения имени) для пользовательского дистрибутива.
phat = mle(ReadmissionTime,'pdf',custpdf,'cdf',custcdf, ... 'Start',0.05,'Censoring',Censored)
phat = 0.1096
Сгенерируйте дважды подвергнутые цензуре данные о выживании и найдите MLEs для встроенного распределения данных. Затем используйте MLEs, чтобы создать объект вероятностного распределения.
Сгенерируйте времена отказа от распределения Бирнбаума-Сондерса.
rng('default') % For reproducibility failuretime = random('BirnbaumSaunders',0.3,1,[100,1]);
Примите, что исследование запускается во время 0.1 и заканчивается во время 0.9. Предположение подразумевает, что времена отказа, меньше чем 0,1 оставляют подвергнутыми цензуре, и времена отказа, больше, чем 0,9, правильные подвергнутый цензуре.
Создайте вектор, в котором каждый элемент указывает на состояние цензуры соответствующего наблюдения в failuretime. Используйте –1, 1, и 0, чтобы указать на лево-подвергнутые цензуре, подвергнутые цензуре правом, и полностью наблюдаемые наблюдения, соответственно.
L = 0.1; U = 0.9; left_censored = (failuretime<L); right_censored = (failuretime>U); c = right_censored - left_censored;
Найдите MLEs для дважды подвергнутых цензуре данных. Укажите информацию цензуры при помощи Censoring аргумент значения имени.
phat = mle(failuretime,'Distribution','BirnbaumSaunders','Censoring',c)
phat = 1×2
0.2632 1.3040
Создайте объект вероятностного распределения с MLEs при помощи makedist функция.
pd = makedist('BirnbaumSaunders','beta',phat(1),'gamma',phat(2))
pd =
BirnbaumSaundersDistribution
Birnbaum-Saunders distribution
beta = 0.263184
gamma = 1.304
pd BirnbaumSaundersDistribution объект. Можно использовать объектные функции pd вычислять распределение и сгенерировать случайные числа. Отобразите поддерживаемые объектные функции.
methods(pd)
Methods for class prob.BirnbaumSaundersDistribution: cdf iqr negloglik proflik truncate gather mean paramci random var icdf median pdf std
Например, вычислите среднее значение и отклонение распределения при помощи mean и var функции, соответственно.
mean(pd)
ans = 0.4869
var(pd)
ans = 0.3681
Сгенерируйте выборочные данные, которые представляют времена аппаратного сбоя после распределения Weibull.
rng('default') % For reproducibility failureTimes = wblrnd(5,2,[200,1]);
Укажите, что наблюдаемые времена отказа являются значениями, округленными к ближайшей секунде.
observed = round(failureTimes);
observed подвергнутые цензуре интервалом данные. Наблюдение t в observed указывает, что событие имело место после времени t–0.5 и перед временем t+0.5.
Создайте матрицу 2D столбца, которая включает информацию о цензуре.
intervalTimes = [observed-0.5 observed+0.5];
Время отказа должно быть положительным. Найдите значения меньшими, чем eps, и измените их в eps.
intervalTimes(intervalTimes < eps) = eps;
Найдите MLEs для параметров распределения Weibull при помощи intervalTimes.
params = mle(intervalTimes,'Distribution','Weibull')
params = 1×2
5.0067 2.0049
Постройте график результатов.
figure histogram(observed,'Normalization','pdf') hold on x = linspace(0,max(observed)); plot(x,wblpdf(x,params(1),params(2))) legend('Observed Samples','Fitted Distribution') hold off
Сгенерируйте выборки от распределения с конечной поддержкой и найдите MLEs с индивидуально настраиваемыми опциями для итеративного процесса оценки.
Для распределения с областью, которая имеет нулевую плотность вероятности, mle может попробовать некоторые параметры, которые имеют нулевую плотность, заставляя функцию не удаться найти MLEs. Чтобы избежать этой проблемы, можно выключить опцию, которая проверяет на недопустимые значения функции, и задайте границы параметра, когда вы вызываете mle функция.
Сгенерируйте выборочные данные размера 1000 от распределения Weibull с масштабным коэффициентом 1 и сформируйте параметр 1. Переключите выборки путем добавления 10.
rng('default') % For reproducibility data = wblrnd(1,1,[1000,1]) + 10; histogram(data,'Normalization','pdf')
Гистограмма не показывает выборок, меньших, чем 10, указывая, что распределение имеет нулевую вероятность в области, меньшей, чем 10. Это распределение является распределением Weibull с тремя параметрами, которое включает третий параметр для местоположения (см. Распределение Weibull С тремя параметрами).
Задайте функцию плотности вероятности (PDF) для распределения Weibull с тремя параметрами.
custompdf = @(x,a,b,c) wblpdf(x-c,a,b);
Найдите MLEs при помощи mle функция. Задайте Options аргумент значения имени, чтобы выключить опцию, которая проверяет на недопустимые значения функции. Кроме того, задайте границы параметра при помощи LowerBound и UpperBound аргументы name-value. Шкала и параметры формы должны быть положительными, и параметр положения должен быть меньшим, чем минимум выборочных данных.
params = mle(data,'pdf',custompdf,'Start',[5 5 5], ... 'Options',statset('FunValCheck','off'), ... 'LowerBound',[0 0 -Inf],'UpperBound',[Inf Inf min(data)])
params = 1×3
1.0258 1.0618 10.0004
mle функция находит точные оценки для этих трех параметров. Для получения дополнительной информации об определении пользовательских опций для итеративного процесса смотрите пример Распределение Weibull С тремя параметрами.
Когда вы предоставляете функции пользовательского дистрибутива или используете встроенные распределения для лево-подвергнутых цензуре, дважды подвергнутых цензуре, подвергнутых цензуре интервалом, или усеченных наблюдений, mle вычисляет оценки параметра с помощью итеративного алгоритма максимизации. С некоторыми моделями и данными, плохим выбором для начальной точки (Start) может вызвать mle сходиться к локальному оптимуму, который не является глобальной переменной maximizer, или не удаваться сходиться полностью. Даже в случаях, для которых логарифмическая правдоподобность хорошего поведения около глобального максимума, выбор начальной точки часто крайне важен для сходимости алгоритма. В частности, если начальные значения параметров далеки от MLEs, потеря значимости в функциях распределения может привести к бесконечной логарифмической правдоподобности.
mle функция находит MLEs путем минимизации отрицательной функции логарифмической правдоподобности (то есть, максимизации функции логарифмической правдоподобности) или при помощи решения закрытой формы, при наличии. Целевая функция является отрицательным значением логарифма продукта выборочных данных (X) вероятности, учитывая параметры распределения (θ):
Функция вероятности P зависит от информации о цензуре для каждого наблюдения.
Полностью наблюдаемое наблюдение — P (x |θ) = f (x), где f является функцией плотности вероятности (PDF) параметрами θ.
Лево-подвергнутое цензуре наблюдение — P (x |θ) = F (x), где F является кумулятивной функцией распределения (cdf) параметрами θ.
Подвергнутое цензуре правом наблюдение — P (x |θ) = 1 – F (x).
Подвергнутое цензуре интервалом наблюдение между xL и xU — P (x |θ) = F (xU) – F (xL).
Для усеченных данных, mle масштабирует функции распределения так, чтобы все вероятности легли в границах усечения [L, U].
mle функция вычисляет доверительные интервалы pci использование точного метода, когда это доступно, и когда выборочные данные не являются усеченными и не включают лево-подвергнутые цензуре или подвергнутые цензуре интервалом наблюдения. В противном случае функция использует Вальдов метод. Точный метод доступен для этих распределений: бином, дискретная универсальная форма, экспоненциальная, нормальная, логарифмически нормальная, Пуассон, Рэлеевская, и непрерывная универсальная форма.