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