Соответствуйте нулевому усеченному распределению Пуассона

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

Во-первых, сгенерируйте некоторые случайные данные Пуассона.

rng(18,'twister') % For reproducibility
lambda = 1.75;
n = 75;
x1 = poissrnd(lambda,n,1);

Затем удалите все нули из данных, чтобы симулировать усечение.

x1 = x1(x1 > 0);

Проверяйте количество отсчетов в x1 после усечения.

length(x1)

ans = 65

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

histogram(x1,0:1:max(x1)+1)

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

Необходимо задать нулевое усеченное распределение Пуассона его функцией вероятностной меры (pmf). Создайте анонимную функцию, чтобы вычислить вероятность для каждой точки в x1, учитывая значение для среднего параметра распределения Пуассона lambda. pmf для нулевого усеченного распределения Пуассона является Пуассон pmf нормированный так, чтобы он суммировал одному. С нулевым усечением нормализацией является 1–Probability(x1<0).

pf_truncpoiss = @(x1,lambda) poisspdf(x1,lambda)./(1-poisscdf(0,lambda));

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

Найдите разумное грубое первое предположение для параметра lambda. В этом случае используйте демонстрационное среднее значение.

start = mean(x1)

start = 2.2154

Обеспечьте mle с данными, пользовательской функцией pmf, начальным значением параметров и нижней границей параметра. Поскольку средний параметр распределения Пуассона должен быть положительным, также необходимо задать нижнюю границу для lambda. mle функция возвращает оценку наибольшего правдоподобия lambda, и опционально, аппроксимированные 95% доверительных интервалов для параметров.

[lambdaHat,lambdaCI] = mle(x1,'pdf',pf_truncpoiss,'Start',start, ...
    'LowerBound',0)

lambdaHat = 1.8760

lambdaCI = 2×1

    1.4990
    2.2530

Оценка параметра меньше, чем демонстрационное среднее значение. Оценка наибольшего правдоподобия составляет нули, не существующие в данных.

В качестве альтернативы можно задать границы усечения при помощи аргумента значения имени TruncationBounds.

[lambdaHat2,lambdaCI2] = mle(x1,'Distribution','Poisson', ...
    'TruncationBounds',[0 Inf])

lambdaHat2 = 1.8760

lambdaCI2 = 2×1

    1.4990
    2.2530

Можно также вычислить оценку стандартной погрешности для lambda при помощи приближения отклонения большой выборки, возвращенного mlecov.

avar = mlecov(lambdaHat,x1,'pdf',pf_truncpoiss);
stderr = sqrt(avar)

stderr = 0.1923

К визуальной проверке подгонка постройте подходящий pmf против нормированной гистограммы необработанных данных

histogram(x1,'Normalization','pdf')
xgrid = min(x1):max(x1);
pmfgrid = pf_truncpoiss(xgrid,lambdaHat);
hold on
plot(xgrid,pmfgrid,'-')
xlabel('x1')
ylabel('Probability')
legend('Sample Data','Fitted pmf','Location','best')
hold off

Соответствуйте верхне-усеченному нормальному распределению

Текущие данные иногда могут b усеченный. Например, наблюдения, больше, чем некоторое фиксированное значение, не могут быть зарегистрированы из-за ограничений в сборе данных.

В этом случае симулируйте данные из усеченного нормального распределения. Во-первых, сгенерируйте некоторые случайные нормальные данные.

n = 500;
mu = 1;
sigma = 3;
rng('default') % For reproducibility
x2 = normrnd(mu,sigma,n,1);

Затем удалите любые наблюдения, которые падают вне точки усечения xTrunc. Примите тот xTrunc известное значение, которое вы не должны оценивать.

xTrunc = 4;
x2 = x2(x2 < xTrunc);

Проверяйте количество отсчетов в x2 после усечения.

length(x2)

ans = 430

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

histogram(x2)

Соответствуйте симулированным данным пользовательским дистрибутивом, который идентичен нормальному распределению для x2 < xTrunc, но имеет нулевую вероятность выше xTrunc. При помощи пользовательского дистрибутива можно оценить нормальные параметры mu и sigma при составлении недостающего хвоста.

Задайте усеченное нормальное распределение его PDF. Создайте анонимную функцию, чтобы вычислить значение плотности вероятности для каждой точки в x, учитывая значения для параметров mu и sigma. С точкой усечения, зафиксированной и известной, PDF усеченным нормальным распределением является усеченная PDF и затем нормированная так, чтобы это объединялось одному. Нормализация является cdf, оцененным в xTrunc. Для простоты примите что весь x2 значения меньше xTrunc, без проверки.

pdf_truncnorm = @(x2,mu,sigma) ...
    normpdf(x2,mu,sigma)./normcdf(xTrunc,mu,sigma);

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

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

start = [mean(x2),std(x2)]

start = 1×2

    0.1585    2.4125

Обеспечьте mle с данными, пользовательской функцией PDF, начальными значениями параметров и нижними границами параметров. Поскольку sigma должно быть положительным, также необходимо задать более низкие границы параметра. mle возвращает оценки наибольшего правдоподобия mu и sigma как один вектор, а также матрица аппроксимированных 95% доверительных интервалов для этих двух параметров.

[paramEsts,paramCIs] = mle(x2,'pdf',pdf_truncnorm,'Start',start, ...
    'LowerBound',[-Inf 0])

paramEsts = 1×2

    1.1298    3.0884

paramCIs = 2×2

    0.5713    2.7160
    1.6882    3.4607

Оценки mu и sigma больше, чем демонстрационное среднее и стандартное отклонение. Подгонка модели составляет недостающий верхний хвост распределения.

В качестве альтернативы можно задать границы усечения при помощи аргумента значения имени TruncationBounds.

[paramEsts2,paramCIs2] = mle(x2,'Distribution','Normal', ...
    'TruncationBounds',[-Inf xTrunc])

paramEsts2 = 1×2

    1.1297    3.0884

paramCIs2 = 2×2

    0.5713    2.7160
    1.6882    3.4607

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

acov = mlecov(paramEsts,x2,'pdf',pdf_truncnorm)

acov = 2×2

    0.0812    0.0402
    0.0402    0.0361

stderr = sqrt(diag(acov))

stderr = 2×1

    0.2849
    0.1900

Чтобы визуально проверять подгонку, постройте подходящую PDF против нормированной гистограммы необработанных данных.

histogram(x2,'Normalization','pdf')
xgrid = linspace(min(x2),max(x2));
pdfgrid = pdf_truncnorm(xgrid,paramEsts(1),paramEsts(2));
hold on
plot(xgrid,pdfgrid,'-')
xlabel('x2')
ylabel('Probability Density')
legend('Sample Data','Fitted pdf','Location','best')
hold off

Подходящая смесь двух нормальных распределений

Некоторые наборы данных показывают bimodality, или даже мультимодальность, и подбор кривой стандартному распределению к таким данным является обычно не соответствующим. Однако смесь простых одномодовых распределений может часто моделировать такие данные очень хорошо.

В этом случае соответствуйте смеси двух нормальных распределений к симулированным данным. Рассмотрите симулированные данные со следующим конструктивным определением:

Сгенерируйте набор данных от смеси t распределений Студента вместо того, чтобы использовать ту же модель, которую вы подбираете. При помощи различных распределений, похожих на метод, используемый при использовании метода Монте-Карло, можно протестировать, насколько устойчивый подходящий метод к отклонениям от предположений о модели, являющейся подходящим.

rng(10) % For reproducibility
x3 = [trnd(20,1,50) trnd(4,1,100)+3];
histogram(x3)

Задайте модель, чтобы соответствовать путем создания анонимной функции, которая вычисляет плотность вероятности. PDF для смеси двух нормальных распределений является взвешенной суммой pdfs двух нормальных компонентов, взвешенных вероятностью смеси. Анонимная функция берет шесть входных параметров: вектор из данных, в которых можно оценить PDF и пять параметров распределения. Каждый компонент имеет параметры для своего среднего и стандартного отклонения.

pdf_normmixture = @(x3,p,mu1,mu2,sigma1,sigma2) ...
    p*normpdf(x3,mu1,sigma1) + (1-p)*normpdf(x3,mu2,sigma2);

Вам также нужно исходное предположение для параметров. Определение начальной точки становится более важным как количество увеличений параметров модели. Здесь, начните с равной смеси (p = 0.5) нормальных распределений, сосредоточенных в двух квартилях данных, с равными стандартными отклонениями. Начальное значение для стандартного отклонения прибывает из формулы для отклонения смеси в терминах среднего значения и отклонения каждого компонента.

pStart = .5;
muStart = quantile(x3,[.25 .75])

muStart = 1×2

    0.3351    3.3046

sigmaStart = sqrt(var(x3) - .25*diff(muStart).^2)

sigmaStart = 1.1602

start = [pStart muStart sigmaStart sigmaStart];

Задайте границы нуля и один для смесительной вероятности и нижних границ нуля для стандартных отклонений. Установите остающиеся элементы векторов границ к +Inf или –Inf, не указать ни на какие ограничения.

lb = [0 -Inf -Inf 0 0];
ub = [1 Inf Inf Inf Inf];
paramEsts = mle(x3,'pdf',pdf_normmixture,'Start',start, ...
    'LowerBound',lb,'UpperBound',ub)

Warning: Maximum likelihood estimation did not converge.  Iteration limit exceeded.

paramEsts = 1×5

    0.3273   -0.2263    2.9914    0.9067    1.2059

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

statset('mlecustom')

ans = struct with fields:
          Display: 'off'
      MaxFunEvals: 400
          MaxIter: 200
           TolBnd: 1.0000e-06
           TolFun: 1.0000e-06
       TolTypeFun: []
             TolX: 1.0000e-06
         TolTypeX: []
          GradObj: 'off'
         Jacobian: []
        DerivStep: 6.0555e-06
      FunValCheck: 'on'
           Robust: []
     RobustWgtFun: []
           WgtFun: []
             Tune: []
      UseParallel: []
    UseSubstreams: []
          Streams: {}
        OutputFcn: []

Максимальное количество по умолчанию итераций для пользовательских дистрибутивов 200. Замените значение по умолчанию, чтобы увеличить число итераций, с помощью структуры опций, созданной с statset функция. Кроме того, увеличьте максимальные вычисления функции.

options = statset('MaxIter',300,'MaxFunEvals',600);
paramEsts = mle(x3,'pdf',pdf_normmixture,'Start',start, ...
    'LowerBound',lb,'UpperBound',ub,'Options',options)

paramEsts = 1×5

    0.3273   -0.2263    2.9914    0.9067    1.2059

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

Чтобы визуально проверять подгонку, постройте подходящую плотность против гистограммы вероятности необработанных данных.

histogram(x3,'Normalization','pdf')
hold on
xgrid = linspace(1.1*min(x3),1.1*max(x3),200);
pdfgrid = pdf_normmixture(xgrid, ...
    paramEsts(1),paramEsts(2),paramEsts(3),paramEsts(4),paramEsts(5));
plot(xgrid,pdfgrid,'-')
hold off
xlabel('x3')
ylabel('Probability Density')
legend('Sample Data','Fitted pdf','Location','best')

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

Mdl = fitgmdist(x3',2)

Mdl = 

Gaussian mixture distribution with 2 components in 1 dimensions
Component 1:
Mixing proportion: 0.329180
Mean:   -0.2200

Component 2:
Mixing proportion: 0.670820
Mean:    2.9975

Mdl.Sigma

ans = 
ans(:,:,1) =

    0.8274


ans(:,:,2) =

    1.4437

Подходящее взвешенное нормальное распределение к данным с неравной точностью

Примите, что у вас есть 10 точек данных, где каждая точка является на самом деле средним значением где угодно от одного до восьми наблюдений. Исходные наблюдения не доступны, но количество наблюдений для каждой точки данных известно. Точность каждой точки зависит от своего соответствующего количества наблюдений. Необходимо оценить среднее и стандартное отклонение необработанных данных.

x4 = [0.25 -1.24 1.38 1.39 -1.43 2.79 3.52 0.92 1.44 1.26]';
m = [8 2 1 3 8 4 2 5 2 4]';

Отклонение каждой точки данных обратно пропорционально своему соответствующему количеству наблюдений, так используйте 1/m взвешивать отклонение каждой точки данных в подгонке наибольшего правдоподобия.

w = 1./m;

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

  c .* exp(-0.5 .* z.^2),

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

Логарифм функции PDF должен вычислить логарифм плотности вероятности для каждой точки в x, учитывая параметры нормального распределения mu и sigma. Этому также нужно с учетом различных весов отклонения. Задайте функцию с именем helper_logpdf_wn1 в отдельном файле helper_logpdf_wn1.m.

function logy = helper_logpdf_wn1(x,m,mu,sigma)
%HELPER_LOGPDF_WN1 Logarithm of pdf for a weight normal distribution
% This function supports only the example Fit Custom Distributions 
% (customdist1demo.mlx) and might change in a future release.
v = sigma.^2 ./ m;
logy = -(x-mu).^2 ./ (2.*v) - .5.*log(2.*pi.*v);
end

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

start = [mean(x4),std(x4)]

start = 1×2

    1.0280    1.5490

Поскольку sigma должно быть положительным, необходимо задать более низкие границы параметра.

[paramEsts1,paramCIs1] = mle(x4,'logpdf', ...
    @(x,mu,sigma)helper_logpdf_wn1(x,m,mu,sigma), ...
    'Start',start,'LowerBound',[-Inf,0])

paramEsts1 = 1×2

    0.6244    2.8823

paramCIs1 = 2×2

   -0.2802    1.6191
    1.5290    4.1456

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

Подходящее нормальное распределение Используя преобразование параметра

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

Во-первых, задайте новый журнал функция с именем PDF helper_logpdf_wn2 это использует преобразованный параметр для sigma.

function logy = helper_logpdf_wn2(x,m,mu,logsigma)
%HELPER_LOGPDF_WN2 Logarithm of pdf for a weight normal distribution with 
% log(sigma) parameterization
% This function supports only the example Fit Custom Distributions 
% (customdist1demo.mlx) and might change in a future release.
v = exp(logsigma).^2 ./ m;
logy = -(x-mu).^2 ./ (2.*v) - .5.*log(2.*pi.*v);
end

Используйте ту же начальную точку, преобразованную для новой параметризации для sigma, то есть, журнал демонстрационного стандартного отклонения.

start = [mean(x4),log(std(x4))]

start = 1×2

    1.0280    0.4376

Поскольку sigma может быть любое положительное значение, log(sigma) неограниченно, и вы не должны задавать нижние или верхние границы.

[paramEsts2,paramCIs2] = mle(x4,'logpdf', ...
    @(x,mu,sigma)helper_logpdf_wn2(x,m,mu,sigma), ...
    'Start',start)

paramEsts2 = 1×2

    0.6244    1.0586

paramCIs2 = 2×2

   -0.2802    0.6203
    1.5290    1.4969

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

sigmaHat = exp(paramEsts2(2))

sigmaHat = 2.8823

sigmaCI = exp(paramCIs2(:,2))

sigmaCI = 2×1

    1.8596
    4.4677

Оценки для обоих mu и sigma эквивалентны в первой подгонке, потому что оценки наибольшего правдоподобия являются инвариантными к параметризации. Доверительный интервал для sigma немного отличается от paramCIs1(:,2).