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

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

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