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

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

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

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.

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

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, но использование log PDF несколько более естественно, потому что обычный PDF имеет форму

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

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

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

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

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

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

start = 1×2

    1.0280    1.5490

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

paramEsts = wgtnormfit

paramEsts = 1×2

    0.6244    2.8823

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

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

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

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

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