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

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

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

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

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

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

ans = 65

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

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

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

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

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

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

start = mean(x)

start = 2.2154

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

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

lambdaHat = 1.8760

lambdaCI = 2×1

    1.4990
    2.2530

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

Мы также можем вычислить стандартную оценку ошибки для лямбды, используя отклонение с большой выборкой приближения возвращенную 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 и sigma с учетом «отсутствующего хвоста».

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

Как и в предыдущем примере, мы зададим усеченное нормальное распределение по его PDF и создадим функцию для вычисления плотности вероятностей для каждой точки x, заданные значения для параметров mu и sigma. С фиксированной и известной точкой усечения 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 и sigma как одного вектора, а также матрицу приблизительных 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 и sigma довольно немного больше, чем среднее значение выборки и стандартное отклонение. Это связано с тем, что подгонка модели рассчитала «отсутствующий» верхний конец распределения.

Мы можем вычислить приблизительную ковариационную матрицу для оценок параметров, используя 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

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

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

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

  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.

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

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

Как и в предыдущих примерах, мы зададим модель для подгонки, создав функцию, которая вычисляет плотность вероятностей. PDF для смеси двух нормалей является всего лишь взвешенной суммой PDF двух нормальных компонентов, взвешенной по вероятности смеси. Этот 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/m, чтобы взвесить отклонение каждой точки данных в максимальной вероятностной подгонке.

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 непосредственно.

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

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

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 in как явного входного параметра.

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

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

start = 1×2

    1.0280    1.5490

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

paramEsts = wgtnormfit

paramEsts = 1×2

    0.6244    2.8823

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

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

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

Новая функция log PDF создается как вложенная функция в функции wgtnormfit2.m. Как и в первой подгонке, этот файл содержит инициализацию данных, вложенную функцию для журнала PDF в взвешенной модели normal и вызов 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 использует ту же начальную точку, преобразованную в новую параметризацию, т.е. журнал стандартного отклонения выборки.

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 и sigma те же, что и в первой подгонке, потому что максимальные оценки правдоподобия инвариантны параметризации.

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