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

Этот пример показывает, как оценить кумулятивную функцию распределения (CDF) от данных непараметрическим или полупараметрическим способом. Это также иллюстрирует метод инверсии для генерации случайных чисел от предполагаемого CDF. Statistics and Machine Learning Toolbox™ включает больше чем две дюжины функций генератора случайных чисел для параметрических одномерных распределений вероятностей. Эти функции позволяют вам генерировать случайные входные параметры для большого разнообразия симуляций, однако, существуют ситуации, где необходимо сгенерировать случайные значения, чтобы моделировать данные, которые не описаны простым параметрическим семейством.

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

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

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

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

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

rng(19,'twister');
n = 25;
x = evrnd(3,1,n,1); x = sort(x);
hist(x,-2:.5:4.5);
xlabel('x'); ylabel('Frequency');

Кусочная линейная непараметрическая оценка CDF

Функция ecdf обеспечивает простой способ вычислить и построить "ступенчатый" эмпирический CDF для данных. В самых простых случаях эта оценка делает дискретные скачки 1/n в каждой точке данных.

[Fi,xi] = ecdf(x);
stairs(xi,Fi,'r');
xlim([-2.5 5]); xlabel('x'); ylabel('F(x)');

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

Просто изменить эмпирический CDF, чтобы обратиться к этому проблемы. Вместо того, чтобы брать дискретные скачки 1/n в каждой точке данных, задайте функцию, которая является кусочна линейный с точками останова в средних точках тех скачков. Высота в каждой из точек данных затем [1/2n, 3/2n..., (n-1/2)/n], вместо [(1/n), (2/n)..., 1]. Используйте вывод ecdf, чтобы вычислить те точки останова, и затем "соединяют точки", чтобы задать кусочную линейную функцию.

xj = xi(2:end);
Fj = (Fi(1:end-1)+Fi(2:end))/2;
hold on
plot(xj,Fj,'b.', xj,Fj,'b-');
hold off
legend({'ECDF' 'Breakpoints' 'Piecewise Linear Estimate'},'location','NW');

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

Поскольку самая маленькая точка данных соответствует высоте 1/2n и самому большому к 1-1/2n, первые и последние линейные сегменты должны быть расширены вне данных, чтобы сделать функциональную досягаемость 0 и 1.

xj = [xj(1)-Fj(1)*(xj(2)-xj(1))/((Fj(2)-Fj(1)));
      xj;
      xj(n)+(1-Fj(n))*((xj(n)-xj(n-1))/(Fj(n)-Fj(n-1)))];
Fj = [0; Fj; 1];
hold on
plot(xj,Fj,'b-','HandleVisibility','off');
hold off

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

F = @(y) interp1(xj,Fj,y,'linear','extrap');
y = linspace(-1,3.75,10);
plot(xj,Fj,'b-',y,F(y),'ko');
xlim([-2.5 5]); xlabel('x'); ylabel('F(x)');

Кусочная линейная непараметрическая обратная оценка CDF

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

stairs(Fi,[xi(2:end); xi(end)],'r');
hold on
plot(Fj,xj,'b-');
hold off
ylim([-2.5 5]); ylabel('x'); xlabel('F(x)');
legend({'ECDF' 'Piecewise Linear Estimate'},'location','NW');

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

Finv = @(u) interp1(Fj,xj,u,'linear','extrap');
u = rand(10000,1);
hist(Finv(u),-2:.25:4.5);
xlabel('x'); ylabel('Frequency');

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

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

Средства оценки ядра для CDF и обратного CDF

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

F = ksdensity(x, x, 'function','cdf', 'width',.35);
stairs(xi,Fi,'r');
hold on
plot(x,F,'b.');
hold off
xlim([-2.5 5]); xlabel('x'); ylabel('F(x)');
legend({'ECDF' 'Kernel Estimates'},'location','NW');

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

y = linspace(-2.5,5,1000);
Fy = ksdensity(x, y, 'function','cdf', 'width',.35);
stairs(xi,Fi,'r');
hold on
plot(y,Fy,'b-');
hold off
legend({'ECDF' 'Kernel Estimate'},'location','NW');

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

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

stairs(Fi,[xi(2:end); xi(end)],'r');
hold on
plot(Fy,y,'b-');
hold off
ylim([-2.5 5]); ylabel('x'); xlabel('F(x)');
legend({'ECDF' 'Kernel Estimate'},'location','NW');

Чтобы преобразовать универсальные случайные значения назад к шкале исходных данных, интерполируйте использование сетки оценок CDF.

Finv = @(u) interp1(Fy,y,u,'linear','extrap');
hist(Finv(u),-2.5:.25:5);
xlabel('x'); ylabel('Frequency');

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

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

r = ksdensity(x, u, 'function','icdf', 'width',.35);

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

r = datasample(x,100000,'replace',true) + normrnd(0,.35,100000,1);

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

binwidth = .25;
edges = -2.5:binwidth:6;
ctrs = edges(1:end-1) + binwidth./2;
counts = histc(r,edges); counts = counts(1:end-1);
bar(ctrs,counts./(sum(counts).*binwidth),1,'FaceColor',[.9 .9 .9]);
hold on
xgrid = edges(1):.1:edges(end);
fgrid = ksdensity(x, xgrid, 'function','pdf', 'width',.3);
plot(xgrid,fgrid,'k-');
hold off
xlabel('x'); ylabel('f(x)');

Полупараметрическая оценка CDF

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

Подбор кривой полупараметрической модели к вашим данным с помощью функции paretotails позволяет лучший из и из непараметрических и из параметрических миров. В "центре" распределения модель использует кусочную линейную непараметрическую оценку для CDF. В каждом хвосте это использует обобщенное распределение Парето. Обобщенный Парето часто используется в качестве модели для хвоста (хвостов) набора данных, и в то время как это достаточно гибко, чтобы соответствовать большому разнообразию хвостов распределения, это достаточно ограничивается так, чтобы это потребовало, чтобы немного данных обеспечили вероятную и сглаженную подгонку к данным о хвосте.

Например, вы можете задать "центр" данных как средние 50% и указать, что переходы между непараметрической оценкой и подгонками Парето происходят в.25 и.75 квантилях ваших данных. Чтобы оценить CDF полупараметрической образцовой подгонки, используйте метод cdf подгонки.

semipFit = paretotails(x,.25,.75);
Warning: Problem fitting generalized Pareto distribution to upper tail. Maximum
likelihood has converged to a boundary point of the parameter space. 

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

[p,q] = boundary(semipFit);
y = linspace(-4,6,1000);
Fy = cdf(semipFit,y);
plot(y,Fy,'b-', q,p,'k+');
xlim([-6.5 5]); xlabel('x'); ylabel('F(x)');
legend({'Semi-parametric Estimate' 'Segment Boundaries'},'location','NW');

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

r = icdf(semipFit,u);
hist(r,-6.5:.25:5);
xlabel('x'); ylabel('Frequency');

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

Заключения

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