Ядерная оценка плотности

Материал из Википедии — свободной энциклопедии
Перейти к навигации Перейти к поиску
Ядерная оценка плотности 100 нормально распределённых случайных чисел с использованием различных сглаживающих полос пропускания.

Ядерная оценка плотности (ЯОП, англ. Kernel Density Estimation, KDE) — это непараметрический способ оценки[англ.] плотности случайной величины. Ядерная оценка плотности является задачей сглаживания данных, когда делается заключение о совокупности, основываясь на конечных выборках данных. В некоторых областях, таких как обработка сигналов и математическая экономика, метод называется также методом окна Парзена-Розенблатта. Как считается, Эммануэль Парзен и Мюррей Розенблатт независимо создали метод в существующем виде[1][2].

Определение

[править | править код]

Пусть является одномерной выборкой независимых одинаково распределённых величин, извлечённых из некоторого распределения с неизвестной плотностью ƒ. Наша задача заключается в оценке формы функции ƒ. Её ядерный оценщик плотности равен

где K является ядром, то есть неотрицательной функцией, а h > 0 является сглаживающим параметром, называемым шириной полосы. Ядро с индексом h называется взвешенным ядром и определяется как . Интуитивно стараются выбрать h как можно меньше, насколько данные это позволяют, однако всегда существует выбор между смещением оценщика и его дисперсией. Выбор полосы пропускания обсуждается более подробно ниже.

Существует ряд наиболее часто используемых ядерных функций: однородная, треугольная, бивзвешенная, тривзвешенная, Епанечникова, нормальная и другие. Ядро Епанечникова оптимально в смысле среднеквадратичной ошибки[3], хотя потеря эффективности для ядер, перечисленных до него, мала[4]. Вследствие удобных математических свойств часто используется нормальное ядро, среднее которого , где является стандартной нормальной функцией плотности.

Построение ядерной оценки плотности находит интерпретацию в областях за пределами оценки плотности[5]. Например, в термодинамике это эквивалентно количеству теплоты, получающейся, когда ядра оператора теплопроводности[англ.] (фундаментальные решения уравнения теплопроводности) помещаются в каждой точке данных xi. Похожие методы используются для построения дискретных операторов Лапласа в точках облака для обучения на основе многообразий[англ.].

Ядерные оценки плотности тесно связаны с гистограммами, но могут быть наделены свойствами, такими как гладкость или непрерывность, путём выбора подходящего ядра. Чтобы это увидеть, сравним построение гистограммы и ядерной оценки плотности на этих 6 точках:

1 2 3 4 5 6
-2,1 -1,3 -0,4 1,9 5,1 6,2

Для гистограммы горизонтальная ось разделена на подинтервалы, которые покрывают область данных. В этом случае мы имеем 6 отрезков, каждая длины 2. Когда точка данных попадает внутрь отрезка, мы помещаем прямоугольник высоты 1/12. Если в отрезок попадает более одной точки, мы размещаем прямоугольники друг над другом.

Для ядерной оценки плотности мы помещаем нормальное ядро с дисперсией 2,25 (показаны красными пунктирными линиями) для каждой точки данных xi. Ядра суммируются, давая ядерную оценку плотности (сплошная синяя кривая). Гладкость ядерной оценки плотности очевидна при сравнении с дискретностью гистограммы, так как ядерные оценки плотности сходятся быстрее к истинной лежащей в основе плотности для непрерывных случайных величин[6].

Сравнение гистограммы (слева) и ядерной оценки плотности (справа), построенных из тех же самых данных. 6 индивидуальных ядер показаны красными пунктирными линиями, ядерная оценка плотности показана синей кривой. Точки данных показаны чёрточками на ленточной диаграмме по горизонтальной оси.
Сравнение гистограммы (слева) и ядерной оценки плотности (справа), построенных из тех же самых данных. 6 индивидуальных ядер показаны красными пунктирными линиями, ядерная оценка плотности показана синей кривой. Точки данных показаны чёрточками на ленточной диаграмме по горизонтальной оси.

Выбор полосы пропускания

[править | править код]
Ядерная оценка плотности (ЯОП, англ. Kernel density estimate, KDE) с различными полосами пропускания случайной выборки 100 точек из стандартного нормального распределения.
Серая кривая: истинная плотность (стандартное нормальное распределение).
Красная кривая: KDE с h=0,05.
Чёрная кривая: KDE с h=0,337.
Зелёная кривая: KDE с h=2.

Полоса пропускания ядра является свободным параметром, который оказывает сильное влияние на результат оценки. Чтобы показать этот эффект, мы возьмём псевдослучайную выборку из обычного нормального распределения (отражена синими чёрточками на ленточной диаграмме[англ.] на горизонтальной оси). Серая кривая представляет истинную плотность (нормальная плотность со средним 0 и дисперсией 1). Для сравнения, красная кривая недостаточно сглажена, поскольку она содержит слишком много случайных выбросов, возникающих при использовании полосы пропускания h=0,05, которая слишком мала. Зелёная кривая чрезмерно сглажена, поскольку используемая полоса пропускания h=2 существенно скрывает структуру. Чёрная кривая с полосой пропускания h=0,337 считается оптимально сглаженной, поскольку её оценка плотности близка к истинной плотности.

Наиболее употребительным критерием оптимальности для выбора этого параметра является ожидаемая функция потерь L2, также называемая средним накопленным квадратом ошибки[англ.] (англ. Mean Integrated Squared Error, MISE):

При слабых предположениях о функциях ƒ и K (ƒ в общем случае является неизвестной вещественной функцией плотности)[1][2], MISE (h)=AMISE(h) + o(1/(nh) + h4), где o является «o» малым. AMISE обозначает «Asymptotic MISE» (асимптотическую MISE), которая состоит из двух ведущих членов

где для функции g, , а ƒ'' является второй производной от ƒ. Для нахождения значения hAMISE, где достигается минимум AMISE, необходимо продифференцировать предыдущее выражение для AMISE по h и получить решение из следующего алгебраического уравнения[7]:

или

Формулы для вычисления AMISE и hAMISE не могут быть использованы прямо, поскольку они включают в себя неизвестную функцию плотности ƒ или её вторую производную ƒ'', так что было разработано большое число автоматических, основанных на данных, методов для выбора полосы пропускания. Во многих обзорах сравнивается эффективность этих методов[8][9][10][11][12][13][14] с общим мнением, что подключаемые выборочные функции[5][15] и функции перекрёстной валидации[16][17][18] наиболее полезны над широкой областью множеств данных.

Подстановка любой полосы пропускания h, которая имеет тот же асимптотический порядок n−1/5 как hAMISE в AMISE даёт , где O — «O» большое. Можно показать, что при слабых предположениях не может существовать непараметрического оценщика, который сходится быстрее, чем ядерный оценщик[19]. Заметим, что скорость n−4/5 меньше, чем типичная скорость сходимости n−1 параметрических методов.

Если полоса пропускания не фиксирована и может меняться в зависимости от места либо величины оценки (balloon оценщик) или величины выборки (поточечный оценщик), получается мощный метод, называемый методом адаптивной ядерной оценки плотности[англ.].

Выбор полосы пропускания для ядерной оценки плотности с медленно убывающим «хвостом» является относительно трудной задачей[20].

Эмпирическое правило для выбора полосы пропускания

[править | править код]

Если используются базовые функции Гаусса для аппроксимации одномерных данных и оцениваемая низлежащая плотность является гауссовой, оптимальный выбор для h (то есть полоса пропускания, минимизирующая средний накопленный квадрат ошибки[англ.]) равна[21]

где является среднеквадратическим отклонением выборки. Аппроксимация называется аппроксимацией нормального распределения, гауссовым распределением или эмпирическим правилом Сильвермана (1986). Хотя это эмпирическое правило вычислительно легко применять, его следует применять с осторожностью, так как оно даёт сильно неточные оценки, когда плотность не близка к нормальной. Например, рассмотрим оценку бимодальной гауссовой смеси:

из выборки с 200 точками. Рисунок справа внизу показывает истинную плотность и две ядерные оценки плотности — одна использует эмпирическое правило выбора полосы, а другая использует выбор полосы, основанный на решении уравнения[5][15]. Оценка на основе эмпирического правила чрезмерно сглажена. Matlab скрипт для примера использует kde.m и дан ниже.

Сравнение между эмпирическим правилом и решением уравнения
Сравнение между эмпирическим правилом и решением уравнения.
% Data
randn('seed',1)                            % Used for reproducibility
data=[randn(100,1)-10; randn(100,1)+10]; % Смесь двух нормальных распределений
% True
phi=@(x) exp(-.5*x.^2)/sqrt(2*pi);       % Нормальная плотность
tpdf=@(x) phi(x+10)/2+phi(x-10)/2;       % Истинная плотность
% Kernel
h=std(data)*(4/3/numel(data))^(1/5);     % Полоса пропускания по эмпирическому правилу Сильвермана
kernel=@(x) mean(phi((x-data)/h)/h);     % Ядерная плотность
kpdf=@(x) arrayfun(kernel,x);            % Поэлементное применение
% Plot
figure(2), clf, hold on
x=linspace(-25,+25,1000);                % Линейная плотность
plot(x,tpdf(x))                            % График истинной плотности
plot(x,kpdf(x))                            % График ядерной плотности с эмпирическим правилом
kde(data)                                  % График ядерной плотности с решением уравнения для вычисления полосы

Связь с характеристической функцией оценщика плотности

[править | править код]

Если дана выборка , естественно оценить характеристическую функцию как

Зная характеристическую функцию, можно найти соответствующую плотность вероятности через формулы преобразования Фурье. Имеется одна трудность применения этой формулы обращения, заключающаяся в том, что это приводит к расходящемуся интегралу, поскольку оценка недостоверна для больших t. Чтобы избежать этой проблемы, оценщик умножается на демпфирующую функцию , которая равна 1 в начале координат, а затем падает до 0 на бесконечности. «Параметр полосы пропускания» h контролирует, насколько мы пытаемся ограничить изменение функции . В частности, когда h мало, будет примерно равно единице для больших t, что означает, что остаётся практически неизменной в наиболее важной области t.

Наиболее употребимым способом выбора функции является либо однородной функцией , которая эффективно означает усечение интервала интегрирования в формуле инвертирования до [−1/h, 1/h], или гауссовой функцией . Когда функция выбрана, может быть применена формула инвертирования и оценщиком плотности будет

где K является преобразованием Фурье демпфирующей функции . Тогда ядерный оценщик плотности совпадает с характеристической функцией оценщика плотности.

Статистические реализации

[править | править код]

Неполный список программного обеспечения, реализующих ядерные оценщики плотности:

  • В пакете Analytica[англ.] релиз 4.4 опция Smoothing для функции плотности вероятности использует KDE, и для выражений возможность доступна как встроенная Pdf функция.
  • На языках C/C++ FIGTree является библиотекой, которая может быть использована для вычисления ядерной оценки плотности с помощью нормальных ядер. Доступен MATLAB-интерфейс.
  • На языке C++ libagf является библиотекой для адаптивной ядерной оценки плотности[англ.].
  • В программе CrimeStat[англ.] ядерная оценка плотности реализована с пятью различными ядерными функциями — нормальная, однородная, четвёртого порядка, отрицательно экспоненциальная и треугольная. Доступны процедуры одно- и двуядерной оценки плотности. Ядерная оценка плотности используется также в интерполяционной процедуре Head Bang, в оценке двухмерной функции плотности Journey-to-crime, и оценке трёхмерной байесовой оценке Journey-to-crime.
  • Во фреймворке ELKI[англ.] ядерные функции плотности можно найти в пакете de.lmu.ifi.dbs.elki.math.statistics.kernelfunctions
  • В продуктах компании ESRI ядерное отображение плотности находится в пакете средств Spatial Analyst и использует ядро четвёртого порядка (бивзвешенное).
  • Для программы Excel компания «the Royal Society of Chemistry» создала дополнение для выполнения ядерной оценки плотности, базирующееся на Analytical Methods Committee Technical Brief 4 (Техническая сводка 4 Комитета Аналитических методов).
  • В gnuplot ядерная оценка плотности реализована опцией smooth kdensity, файл данных может содержать вес и полосу пропускания для каждой точки или полоса пропускания может быть установлена автоматически[22] согласно «эмпирическому правилу Сильвермана» (см. выше).
  • В языке Haskell ядерная плотность реализована в пакете statistics.
  • В системе IGOR Pro[англ.] ядерная оценка плотности реализована в виде операции StatsKDE (добавлена в версии Igor Pro 7.00). Полоса пропускания может быть указана или оценена средними Сильвермана, Скотта или Боуманна и Аззалини. Типы ядер: Епанечникова, бивзвешенные, тривзвешенные, треугольные, гауссовы и прямоугольные.
  • На языке Java, пакет Weka предоставляет weka.estimators.KernelEstimator среди прочего другого.
  • На языке JavaScript пакет визуализации D3.js[англ.] содержит пакет KDE в пакете science.stats.
  • В пакете JMP[англ.] может быть использовано средство «Distribution platform» для создания одномерной ядерной оценки плотности, а «Fit Y by X platform» может быть использовано для создания двумерной ядерной оценки плотности.
  • На языке Julia ядерная оценка плотности реализована в пакете KernelDensity.jl.
  • В программе MATLAB ядерная оценка плотности реализована через функцию ksdensity (Statistics Toolbox). В релизе 2018 года MATLAB, могут быть заданы как полоса пропускания, так и ядерный сглаживатель, включая и другие опции, такие как указание пределов ядерной плотности. Альтернативно, бесплатный пакет для MATLAB, который реализует автоматический выбор полосы пропускания[5], доступен на странице «MATLAB Central File Exchange» для
  • В системе Mathematica численная оценка ядерного распределения реализована в виде функции SmoothKernelDistribution здесь, а символьная оценка реализована с помощью функции KernelMixtureDistribution here и обе реализации осуществляют выбор полосы пропускания из представленных данных.
  • Для пакета Minitab компания «the Royal Society of Chemistry» создала макро для осуществления ядерной оценки плотности на основе их Analytical Methods Committee Technical Brief 4 (Техническая сводка 4 Комитета Аналитических методов).
  • В библиотеке NAG[англ.] ядерная оценка плотности реализована процедурой g10ba (доступной на языке Fortran[24] и языке C[25]).
  • В библиотеке Nuklei методы ядерной плотности на языке C++ фокусируются на дынных из специальной евклидовой группы .
  • В системе Octave ядерная оценка плотности реализована возможностью kernel_density (пакет математической экономики).
  • В пакете Origin 2D ядерный график плотности может быть построен с помощью пользовательского интерфейса пакета, а коды двух функций Ksdensity для 1D и Ks2density для 2D могут быть взяты на языках LabTalk, Python или C.
  • На языке Perl реализацию можно найти в модуле Statistics-KernelEstimati
  • На языке PHP реализацию можно найти в библиотеке MathPHP
  • На языке Python существует множество реализаций: pyqt_fit.kde Module в пакете PyQt-Fit, SciPy (scipy.stats.gaussian_kde и scipy.signal.parzen), Statsmodels (KDEUnivariate и KDEMultivariate) и Scikit-learn (KernelDensity) (см. сравнение[26]). KDEpy поддерживает взвешенные данные, и FFT-реализация на порядок быстрее других реализаций.
  • В языке R это реализовано через density в базовой поставке, через bkde в библиотеке KernSmooth, через ParetoDensityEstimation в библиотеке AdaptGauss (для оценки плотности распределения Парето), через kde в библиотеке ks, через dkden и dbckden в библиотеке evmix, npudens в библиотеке np (численные и категоральные данные), sm.density в библиотеке sm. Для получения реализации функции kde.R, которая не требует установки какого-либо пакета или библиотеки, см. kde.R. Библиотека btb, предназначенная для городского анализа, реализует ядерную оценку плотности через kernel_smoothing.
  • В системе SAS (программа)[англ.] может быть использована процедура proc kde для оценки одномерных и двумерных ядерных плотностей.
  • В пакете Stata[англ.] это реализовано в виде kdensity[27], например histogram x, kdensity. Альтернативно, доступен бесплатный модуль KDENS пакета Stata здесь, который позволяет оценить 1D- или 2D-функции плотности.
  • В Apache Spark вы можете использовать класс KernelDensity() (см. официальную документацию)

Примечания

[править | править код]
  1. 1 2 Rosenblatt, 1956, с. 832.
  2. 1 2 Parzen, 1962, с. 1065.
  3. Epanechnikov, 1969, с. 153–158.
  4. Wand, Jones, 1995.
  5. 1 2 3 4 Botev, Grotowski, Kroese, 2010, с. 2916–2957.
  6. Scott, 1979, с. 605–610.
  7. В. А. Епанечников, “Непараметрическая оценка многомерной плотности вероятности”, Теория вероятн. и ее примен., 14:1 (1969), 156–161; Theory Probab. Appl., 14:1 (1969), 153–158. www.mathnet.ru. Дата обращения: 31 января 2022.
  8. Park, Marron, 1990, с. 66–72.
  9. Park, Turlach, 1992, с. 251–270.
  10. Cao, Cuevas, Manteiga, 1994, с. 153–176.
  11. Jones, Marron, Sheather, 1996, с. 401–407.
  12. Sheather, 1992, с. 225–250, 271–281.
  13. Agarwal, Aluru, 2010, с. 575–597.
  14. Xu, Yan, Xu, 2015, с. 28–37.
  15. 1 2 Sheather, Jones, 1991, с. 683–690.
  16. Rudemo, 1982, с. 65–78.
  17. Bowman, 1984, с. 353–360.
  18. Hall, Marron, Park, 1992, с. 1–20.
  19. Wahba, 1975, с. 15–29.
  20. Buch-Larsen, 2005, с. 503–518.
  21. Silverman, 1986, с. 48.
  22. Janert, 2009, с. section 13.2.2.
  23. Horová, Koláček, Zelinka, 2012.
  24. The Numerical Algorithms Group NAG Library Routine Document: nagf_smooth_kerndens_gauss (g10baf). NAG Library Manual, Mark 23. Дата обращения: 16 февраля 2012. Архивировано 3 апреля 2015 года.
  25. The Numerical Algorithms Group NAG Library Routine Document: nag_kernel_density_estim (g10bac). NAG Library Manual, Mark 9. Дата обращения: 16 февраля 2012. Архивировано из оригинала 24 ноября 2011 года.
  26. Vanderplas, Jake Kernel Density Estimation in Python (1 декабря 2013). Дата обращения: 12 марта 2014.
  27. Источник. Дата обращения: 27 октября 2018. Архивировано 29 августа 2017 года.

Литература

[править | править код]
  • Rosenblatt M. Remarks on Some Nonparametric Estimates of a Density Function // The Annals of Mathematical Statistics. — 1956. — Т. 27, вып. 3. — doi:10.1214/aoms/1177728190.
  • Parzen E. On Estimation of a Probability Density Function and Mode // The Annals of Mathematical Statistics. — 1962. — Т. 33, вып. 3. — doi:10.1214/aoms/1177704472. — JSTOR 2237880.
  • Epanechnikov V.A. Non-parametric estimation of a multivariate probability density // Theory of Probability and its Applications. — 1969. — Т. 14. — doi:10.1137/1114019.
  • Wand M.P, Jones M.C. Kernel Smoothing. — London: Chapman & Hall/CRC, 1995. — ISBN 0-412-55270-1.
  • Botev Z.I., Grotowski J.F., Kroese D.P. Kernel density estimation via diffusion // Annals of Statistics. — 2010. — Т. 38, вып. 5. — doi:10.1214/10-AOS799. — arXiv:1011.2602.
  • Scott D. On optimal and data-based histograms // Biometrika. — 1979. — Т. 66, вып. 3. — doi:10.1093/biomet/66.3.605.
  • Park B. U., Marron J. S. Comparison of data-driven bandwidth selectors // Journal of the American Statistical Association. — 1990. — Т. 85, вып. 409. — doi:10.1080/01621459.1990.10475307. — JSTOR 2289526.
  • Park B. U., Turlach B. A. Practical performance of several data driven bandwidth selectors (with discussion) // Computational Statistics. — 1992. — Т. 7. — С. 251–270.
  • Cao R., Cuevas A., Manteiga W. G. A comparative study of several smoothing methods in density estimation // Computational Statistics and Data Analysis. — 1994. — Т. 17, вып. 2. — doi:10.1016/0167-9473(92)00066-Z.
  • Jones M.C., Marron J.S., Sheather S. J. A brief survey of bandwidth selection for density estimation // Journal of the American Statistical Association. — 1996. — Т. 91, вып. 433. — doi:10.2307/2291420. — JSTOR 2291420.
  • Sheather S.J. The performance of six popular bandwidth selection methods on some real data sets (with discussion) // Computational Statistics. — 1992. — Т. 7.
  • Agarwal N., Aluru N.R. A data-driven stochastic collocation approach for uncertainty quantification in MEMS // International Journal for Numerical Methods in Engineering. — 2010. — Т. 83, вып. 5.
  • Xu X., Yan Z., Xu S. Estimating wind speed probability distribution by diffusion-based kernel density method // Electric Power Systems Research. — 2015. — Т. 121. — С. 28–37.
  • Sheather S.J., Jones M.C. A reliable data-based bandwidth selection method for kernel density estimation // Journal of the Royal Statistical Society, Series B. — 1991. — Т. 53, вып. 3. — JSTOR 2345597.
  • Rudemo M. Empirical choice of histograms and kernel density estimators // Scandinavian Journal of Statistics. — 1982. — Т. 9, вып. 2. — JSTOR 4615859.
  • Bowman A.W. An alternative method of cross-validation for the smoothing of density estimates // Biometrika. — 1984. — Т. 71, вып. 2. — doi:10.1093/biomet/71.2.353.
  • Hall P., Marron J.S., Park B.U. Smoothed cross-validation // Probability Theory and Related Fields. — 1992. — Т. 92. — С. 1–20. — doi:10.1007/BF01205233.
  • Wahba G. Optimal convergence properties of variable knot, kernel, and orthogonal series methods for density estimation // Annals of Statistics. — 1975. — Т. 3, вып. 1. — doi:10.1214/aos/1176342997.
  • TINE Buch-Larsen. Kernel density estimation for heavy-tailed distributions using the Champernowne transformation // Statistics. — 2005. — Т. 39, вып. 6. — doi:10.1080/02331880500439782.
  • Silverman B.W. Density Estimation for Statistics and Data Analysis. — London: Chapman & Hall/CRC, 1986. — ISBN 0-412-24620-1.
  • Philipp K Janert. section 13.2.2 Kernel density estimates // Gnuplot in action : understanding data with graphs. — Connecticut, USA: Manning Publications, 2009. — ISBN 978-1-933988-39-9.
  • Horová I., Koláček J., Zelinka J. Kernel Smoothing in MATLAB: Theory and Practice of Kernel Smoothing. — Singapore: World Scientific Publishing, 2012. — ISBN 978-981-4405-48-5.
  • Introduction to kernel density estimation Короткое введение, где объясняется ядерная оценка плотности как усовершенствование гистограмм.
  • Kernel Bandwidth Optimization Свободное онлайн средство, которое генерирует оптимизированную ядерную оценку плотности на ваших данных.
  • Free Online Software (Calculator) вычисляет ядерную оценку плотности для любой выборки для ядер: гауссово, Епанечникова, прямоугольное, треугольное, бивзвешенное, косинусное и optcosine.
  • Kernel Density Estimation Applet Онлайновый интерактивный пример ядерной оценки плотности. Требует версию NET 3.0 или выше.