Ядерная оценка плотности (ЯОП, , KDE) — непараметрический способ случайной величины. Представляет собой задачу сглаживания данных, в результате решения которой на основании конечных выборок данных делается заключение о генеральной совокупности. В некоторых областях, таких как обработка сигналов и математическая экономика, ЯОП называется также методом окна Парзена — Розенблатта. Считается, что и Мюррей Розенблатт независимо создали метод в существующем виде.
Определение Пусть — одномерная выборка независимых одинаково распределённых величин, извлечённых из некоторого распределения с неизвестной плотностью , которую необходимо оценить. Её ядерный оценщик плотности равен
где — ядро, то есть некая неотрицательная функция, — сглаживающий параметр, называемый шириной полосы. Ядро с индексом называется взвешенным ядром и определяется как . С точки зрения интуиции, нужно выбирать настолько маленьким, насколько это позволяют данные, однако на практике всегда приходится искать компромисс между смещением оценки и её дисперсией .
Существует ряд наиболее часто используемых ядерных функций: однородная, треугольная, бивзвешенная, тривзвешенная, Епанечникова, нормальная и другие. Ядро Епанечникова оптимально в смысле среднеквадратичной ошибки, хотя для остальных перечисленных ядер потеря эффективности мала. Вследствие присущих ему удобных математических свойств, часто используется нормальное ядро, для которого , где — функция плотности стандартного нормального распределения.
Построение ядерной оценки плотности находит свои интерпретации в областях за пределами оценки плотности. Например, в термодинамике ядерная оценка эквивалентна количеству теплоты, получающейся, когда (фундаментальные решения уравнения теплопроводности) помещаются в каждой точке данных . Похожие методы используются для построения дискретных операторов Лапласа на точечных облаках для .
Ядерные оценки плотности тесно связаны с гистограммами, но могут быть наделены свойствами, такими как гладкость или непрерывность, путём выбора подходящего ядра. Чтобы убедиться в этом, сравним построение гистограммы и ядерной оценки плотности на следующих 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 (показаны красными пунктирными линиями) для каждой точки данных . Ядра суммируются, давая ядерную оценку плотности (сплошная синяя кривая). Гладкость ядерной оценки плотности очевидна при сравнении с дискретностью гистограммы, так как ядерные оценки плотности сходятся быстрее к истинной плотности, лежащей в основе данных, для непрерывных случайных величин.
Выбор полосы пропускания
Полоса пропускания ядра является свободным параметром, который оказывает сильное влияние на результат оценки. Чтобы показать этот эффект, мы возьмём псевдослучайную выборку из обычного нормального распределения (отражена синими чёрточками на на горизонтальной оси). Серая кривая представляет истинную плотность (нормальная плотность со средним 0 и дисперсией 1). Для сравнения, красная кривая недостаточно сглажена, поскольку она содержит слишком много случайных выбросов, возникающих при использовании полосы пропускания h=0,05, которая слишком мала. Зелёная кривая чрезмерно сглажена, поскольку используемая полоса пропускания h=2 существенно скрывает структуру. Чёрная кривая с полосой пропускания h=0,337 считается оптимально сглаженной, поскольку её оценка плотности близка к истинной плотности.
Наиболее употребительным критерием оптимальности для выбора этого параметра является ожидаемая функция потерь L2, также называемая (, MISE):
При слабых предположениях о функциях ƒ и K (ƒ в общем случае является неизвестной вещественной функцией плотности), MISE (h)=AMISE(h) + o(1/(nh) + h4), где o является «o» малым. AMISE обозначает «Asymptotic MISE» (асимптотическую MISE), которая состоит из двух ведущих членов
где для функции g, , а ƒ'''' является второй производной от ƒ. Для нахождения значения hAMISE, где достигается минимум AMISE, необходимо продифференцировать предыдущее выражение для AMISE по h и получить решение из следующего алгебраического уравнения:
или
Формулы для вычисления AMISE и hAMISE не могут быть использованы прямо, поскольку они включают в себя неизвестную функцию плотности ƒ или её вторую производную ƒ'', так что было разработано большое число автоматических, основанных на данных, методов для выбора полосы пропускания. Во многих обзорах сравнивается эффективность этих методов с общим мнением, что подключаемые выборочные функции и функции перекрёстной валидации наиболее полезны над широкой областью множеств данных.
Подстановка любой полосы пропускания h, которая имеет тот же асимптотический порядок n−1/5 как hAMISE в AMISE даёт , где O — «O» большое. Можно показать, что при слабых предположениях не может существовать непараметрического оценщика, который сходится быстрее, чем ядерный оценщик. Заметим, что скорость n−4/5 меньше, чем типичная скорость сходимости n−1 параметрических методов.
Если полоса пропускания не фиксирована и может меняться в зависимости от места либо величины оценки (balloon оценщик) или величины выборки (поточечный оценщик), получается мощный метод, называемый методом .
Выбор полосы пропускания для ядерной оценки плотности с медленно убывающим «хвостом» является относительно трудной задачей.
Эмпирическое правило для выбора полосы пропускания Если используются базовые функции Гаусса для аппроксимации одномерных данных и оцениваемая низлежащая плотность является гауссовой, оптимальный выбор для h (то есть полоса пропускания, минимизирующая ) равна
где является среднеквадратическим отклонением выборки. Аппроксимация называется аппроксимацией нормального распределения, гауссовым распределением или эмпирическим правилом Сильвермана (1986). Хотя это эмпирическое правило вычислительно легко применять, его следует применять с осторожностью, так как оно даёт сильно неточные оценки, когда плотность не близка к нормальной. Например, рассмотрим оценку бимодальной гауссовой смеси:
из выборки с 200 точками. Рисунок справа внизу показывает истинную плотность и две ядерные оценки плотности — одна использует эмпирическое правило выбора полосы, а другая использует выбор полосы, основанный на решении уравнения. Оценка на основе эмпирического правила чрезмерно сглажена. 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.
Наиболее употребимым способом выбора функции является либо однородной функцией , которая эффективно означает усечение интервала интегрирования в формуле инвертирования до , или гауссовой функцией . Когда функция выбрана, может быть применена формула инвертирования и оценщиком плотности будет
где K является преобразованием Фурье демпфирующей функции . Тогда ядерный оценщик плотности совпадает с характеристической функцией оценщика плотности.
Статистические реализации Неполный список программного обеспечения, реализующих ядерные оценщики плотности: В пакете релиз 4.4 опция Smoothing'' для функции плотности вероятности использует KDE, и для выражений возможность доступна как встроенная Pdf функция. На языках C/C++ FIGTree является библиотекой, которая может быть использована для вычисления ядерной оценки плотности с помощью нормальных ядер. Доступен MATLAB-интерфейс. На языке C++ libagf является библиотекой для . В программе ядерная оценка плотности реализована с пятью различными ядерными функциями — нормальная, однородная, четвёртого порядка, отрицательно экспоненциальная и треугольная. Доступны процедуры одно- и двуядерной оценки плотности. Ядерная оценка плотности используется также в интерполяционной процедуре Head Bang, в оценке двухмерной функции плотности Journey-to-crime, и оценке трёхмерной байесовой оценке Journey-to-crime. Во фреймворке ядерные функции плотности можно найти в пакете 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, файл данных может содержать вес и полосу пропускания для каждой точки или полоса пропускания может быть установлена автоматически согласно «эмпирическому правилу Сильвермана» (см. выше). В языке Haskell ядерная плотность реализована в пакете statistics. В системе ядерная оценка плотности реализована в виде операции StatsKDE (добавлена в версии Igor Pro 7.00). Полоса пропускания может быть указана или оценена средними Сильвермана, Скотта или Боуманна и Аззалини. Типы ядер: Епанечникова, бивзвешенные, тривзвешенные, треугольные, гауссовы и прямоугольные. На языке Java, пакет Weka предоставляет weka.estimators.KernelEstimator среди прочего другого. На языке JavaScript пакет визуализации содержит пакет KDE в пакете science.stats. В пакете JMP может быть использовано средство «Distribution platform» для создания одномерной ядерной оценки плотности, а «Fit Y by X platform» может быть использовано для создания двумерной ядерной оценки плотности. На языке Julia ядерная оценка плотности реализована в пакете KernelDensity.jl. В программе MATLAB ядерная оценка плотности реализована через функцию ksdensity (Statistics Toolbox). В релизе 2018 года MATLAB, могут быть заданы как полоса пропускания, так и ядерный сглаживатель, включая и другие опции, такие как указание пределов ядерной плотности. Альтернативно, бесплатный пакет для MATLAB, который реализует автоматический выбор полосы пропускания, доступен на странице «MATLAB Central File Exchange» для 1-мерных данных 2-мерных данных n-мерных данные Бесплатный набор средств для MATLAB с реализацией ядерной регрессии, ядерной оценки плотности, ядерной оценки функции риска и многих других доступен на этой странице (этот набор является частью книги «Kernel Smoothing in MATLAB»).