Подбор кривой пользовательским дистрибутивам: нулевой усеченный пример Пуассона

Рассчитайте данные часто моделируются с помощью распределения Пуассона, и можно использовать функцию Statistics and Machine Learning Toolbox poissfit подбирать модель Poisson. Однако в некоторых ситуациях, количества, которые являются нулем, не становятся записанными в данных, и так подбор кривой распределению Пуассона не является прямым из-за тех "недостающих нулей". Этот пример покажет, как соответствовать распределению Пуассона к нулевым усеченным данным, с помощью функционального mle.

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

rng(18,'twister');
n = 75;
lambda = 1.75;
x = poissrnd(lambda,n,1);

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

x = x(x > 0);
length(x)

ans = 65

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

hist(x,[0:1:max(x)+1]);

Первый шаг должен задать нулевое усеченное распределение Пуассона своей функцией вероятности (PF). Мы создадим функцию, чтобы вычислить вероятность для каждой точки в x, учитывая значение для среднего lambda параметра распределения Пуассона. PF для нулевого усеченного Пуассона является только обычным PF Пуассона, повторно нормированным так, чтобы это суммировало одному. С нулевым усечением перенормализация только с 1 PR {0}. Самый легкий способ создать функцию для PF состоит в том, чтобы использовать анонимную функцию.

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

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

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

start = mean(x)

start = 2.2154

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

[lambdaHat,lambdaCI] = mle(x, 'pdf',pf_truncpoiss, 'start',start, 'lower',0)

lambdaHat = 1.8760

lambdaCI = 2×1

    1.4990
    2.2530

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

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

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

stderr = 0.1923

Предоставление дополнительных значений к функции распределения: усеченный нормальный пример

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

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

n = 75;
mu = 1;
sigma = 3;
x = normrnd(mu,sigma,n,1);

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

xTrunc = 4;
x = x(x < xTrunc);
length(x)

ans = 64

Вот гистограмма этих симулированных данных. Мы будем соответствовать им распределением, которое идентично нормальному для x <xTrunc, но это имеет нулевую вероятность выше xTrunc. Таким образом мы можем оценить нормальные параметры mu и сигму при объяснении "недостающего хвоста".

hist(x,(-10:.5:4));

Как в предыдущем примере, мы зададим усеченное нормальное распределение его PDF и создадим функцию, чтобы вычислить плотность вероятности для каждой точки в x, учитывая значения для параметров mu и сигмы. С точкой усечения, зафиксированной и известной, PDF усеченным нормальным является только обычная нормальная PDF, усеченная, и затем повторно нормированная так, чтобы это объединялось одному. Перенормализация является только CDF, оцененным в xTrunc. Для простоты мы примем, что все x значения меньше xTrunc без проверки. Мы будем использовать анонимную функцию, чтобы задать PDF.

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

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

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

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

start = 1×2

    0.2735    2.2660

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

[paramEsts,paramCIs] = mle(x, 'pdf',pdf_truncnorm, 'start',start, 'lower',[-Inf 0])

paramEsts = 1×2

    0.9911    2.7800

paramCIs = 2×2

   -0.1605    1.9680
    2.1427    3.5921

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

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

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

acov = 2×2

    0.3452    0.1660
    0.1660    0.1717

stderr = sqrt(diag(acov))

stderr = 2×1

    0.5876
    0.4143

Подбор кривой более сложному распределению: смесь двух нормалей

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

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

  First, flip a biased coin.  If it lands heads, pick a value at random
  from a normal distribution with mean mu_1 and standard deviation
  sigma_1. If the coin lands tails, pick a value at random from a normal
  distribution with mean mu_2 and standard deviation sigma_2.

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

x = [trnd(20,1,50) trnd(4,1,100)+3];
hist(x,-2.25:.5:7.25);

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

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

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

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

muStart = 1×2

    0.5970    3.2456

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

sigmaStart = 1.2153

start = [pStart muStart sigmaStart sigmaStart];

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

lb = [0 -Inf -Inf 0 0];
ub = [1 Inf Inf Inf Inf];
paramEsts = mle(x, 'pdf',pdf_normmixture, 'start',start, ...
                          'lower',lb, 'upper',ub)

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

paramEsts = 1×5

    0.3523    0.0257    3.0489    1.0546    1.0629

Значением по умолчанию для пользовательских дистрибутивов являются 200 итераций.

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: []

Замените то значение по умолчанию, с помощью структуры опций, созданной с statset функция. Также увеличьте предел вычисления функции.

options = statset('MaxIter',300, 'MaxFunEvals',600);
paramEsts = mle(x, 'pdf',pdf_normmixture, 'start',start, ...
                          'lower',lb, 'upper',ub, 'options',options)

paramEsts = 1×5

    0.3523    0.0257    3.0489    1.0546    1.0629

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

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

bins = -2.5:.5:7.5;
h = bar(bins,histc(x,bins)/(length(x)*.5),'histc');
h.FaceColor = [.9 .9 .9];
xgrid = linspace(1.1*min(x),1.1*max(x),200);
pdfgrid = pdf_normmixture(xgrid,paramEsts(1),paramEsts(2),paramEsts(3),paramEsts(4),paramEsts(5));
hold on
plot(xgrid,pdfgrid,'-')
hold off
xlabel('x')
ylabel('Probability Density')

Используя вложенные функции: нормальный пример с неравной точностью

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

Примите, что у нас есть 10 точек данных, где каждый - на самом деле среднее значение где угодно от 1 до 8 наблюдений. Те исходные наблюдения не доступны, но мы знаем, сколькими было для каждой из наших точек данных. Мы должны оценить среднее и стандартное отклонение необработанных данных.

x = [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/м, чтобы взвесить отклонение каждой точки данных в подгонке наибольшего правдоподобия.

w = 1 ./ m

w = 1×10

    0.1250    0.5000    1.0000    0.3333    0.1250    0.2500    0.5000    0.2000    0.5000    0.2500

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

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

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

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

Мы создали отдельный файл для функции под названием wgtnormfit.m. Эта функция содержит инициализацию данных, вложенную функцию для журнала PDF во взвешенной нормальной модели и вызов mle функционируйте, чтобы на самом деле подбирать модель. Поскольку сигма должна быть положительной, мы должны задать более низкие границы параметра. Вызов mle возвращает оценки наибольшего правдоподобия mu и сигмы в одном векторе.

type wgtnormfit.m

function paramEsts = wgtnormfit
%WGTNORMFIT Fitting example for a weighted normal distribution.

%   Copyright 1984-2012 The MathWorks, Inc.


x = [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]';

    function logy = logpdf_wn(x,mu,sigma)
        v = sigma.^2 ./ m;
        logy = -(x-mu).^2 ./ (2.*v) - .5.*log(2.*pi.*v);
    end

paramEsts = mle(x, 'logpdf',@logpdf_wn, 'start',[mean(x),std(x)], 'lower',[-Inf,0]);

end

В wgtnormfit.m, мы передаем mle указатель на вложенную функцию logpdf_wn, использование 'logpdf' параметра. Та вложенная функция относится к количествам наблюдения, m, в расчете взвешенного журнала PDF. Поскольку вектор m задан в его родительской функции, logpdf_wn имеет доступ к нему, и нет никакой потребности волноваться о передаче m в как явный входной параметр.

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

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

start = 1×2

    1.0280    1.5490

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

paramEsts = wgtnormfit

paramEsts = 1×2

    0.6244    2.8823

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

Используя преобразование параметра: нормальный пример (продолжен)

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

Новый журнал функция PDF создается как вложенная функция в функциональном wgtnormfit2.m. Как в первой подгонке, этот файл содержит инициализацию данных, вложенную функцию для журнала PDF во взвешенной нормальной модели и вызов mle функционируйте, чтобы на самом деле подбирать модель. Поскольку сигма может быть любым положительным значением, журнал (сигма) неограничен, и мы больше не должны задавать нижние или верхние границы. Кроме того, вызов mle в этом случае возвращает и оценки параметра и доверительные интервалы.

type wgtnormfit2.m

function [paramEsts,paramCIs] = wgtnormfit2
%WGTNORMFIT2 Fitting example for a weighted normal distribution (log(sigma) parameterization).

%   Copyright 1984-2012 The MathWorks, Inc.


x = [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]';

    function logy = logpdf_wn2(x,mu,logsigma)
        v = exp(logsigma).^2 ./ m;
        logy = -(x-mu).^2 ./ (2.*v) - .5.*log(2.*pi.*v);
    end

[paramEsts,paramCIs] = mle(x, 'logpdf',@logpdf_wn2, 'start',[mean(x),log(std(x))]);

end

Заметьте тот wgtnormfit2.m использует ту же начальную точку, преобразованную к новой параметризации, i.e., журнал демонстрационного стандартного отклонения.

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

start = 1×2

    1.0280    0.4376

[paramEsts,paramCIs] = wgtnormfit2

paramEsts = 1×2

    0.6244    1.0586

paramCIs = 2×2

   -0.2802    0.6203
    1.5290    1.4969

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

muHat = paramEsts(1)

muHat = 0.6244

sigmaHat = exp(paramEsts(2))

sigmaHat = 2.8823

muCI = paramCIs(:,1)

muCI = 2×1

   -0.2802
    1.5290

sigmaCI = exp(paramCIs(:,2))

sigmaCI = 2×1

    1.8596
    4.4677