Оценка многомерной плотности ядра - Multivariate kernel density estimation

Оценка плотности ядра это непараметрический техника для оценка плотности т.е. оценка функции плотности вероятности, что является одним из основных вопросов в статистика. Его можно рассматривать как обобщение гистограмма оценка плотности с улучшенными статистическими свойствами. Помимо гистограмм, к другим типам оценок плотности относятся: параметрический, сплайн, вейвлет и Ряд Фурье. Оценщики плотности ядра были впервые представлены в научной литературе для одномерный данные в 1950-х и 1960-х годах[1][2] и впоследствии получили широкое распространение. Вскоре было признано, что аналогичные оценки для многомерных данных будут важным дополнением к многомерная статистика. На основании исследований, проведенных в 1990-х и 2000-х годах, многомерная оценка плотности ядра достигла уровня зрелости, сопоставимого с его одномерными аналогами.[3]

Мотивация

Возьмем иллюстративный синтетический двумерный набор данных из 50 точек для иллюстрации построения гистограмм. Для этого требуется выбор точки привязки (нижний левый угол сетки гистограммы). Для гистограммы слева мы выбираем (-1,5, -1,5): для гистограммы справа мы сдвигаем точку привязки на 0,125 в обоих направлениях до (-1,625, -1,625). Обе гистограммы имеют ширину бина 0,5, поэтому любые различия связаны только с изменением точки привязки. Цветовая кодировка указывает количество точек данных, которые попадают в ячейку: 0 = белый, 1 = бледно-желтый, 2 = ярко-желтый, 3 = оранжевый, 4 = красный. Левая гистограмма, по-видимому, указывает на то, что верхняя половина имеет более высокую плотность, чем нижняя половина, тогда как обратное верно для правой гистограммы, подтверждая, что гистограммы очень чувствительны к размещению точки привязки.[4]

Оставили. Histogram with anchor point at (−1.5, -1.5). Правильно. Histogram with anchor point at (−1.625, −1.625). Both histograms have a bin width of 0.5, so differences in appearances of the two histograms are due to the placement of the anchor point.
Сравнение двухмерных гистограмм. Оставили. Гистограмма с точкой привязки в (-1,5, -1,5). Правильно. Гистограмма с точкой привязки в (-1,625, -1,625). Обе гистограммы имеют ширину ячейки 0,5, поэтому различия во внешнем виде двух гистограмм связаны с расположением точки привязки.

Одним из возможных решений этой проблемы размещения точек привязки является полное удаление сетки биннинга гистограммы. На левом рисунке ниже ядро ​​(представленное серыми линиями) центрировано в каждой из 50 точек данных выше. Результат суммирования этих ядер показан на правом рисунке, который является оценкой плотности ядра. Наиболее разительное различие между оценками плотности ядра и гистограммами заключается в том, что первые легче интерпретировать, поскольку они не содержат искажений, вызванных сеткой бинирования. Цветные контуры соответствуют наименьшей области, которая содержит соответствующую вероятностную массу: красный = 25%, оранжевый + красный = 50%, желтый + оранжевый + красный = 75%, что указывает на то, что одна центральная область имеет самую высокую плотность.

Оставили. Individual kernels. Правильно. Kernel density estimate.
Построение 2D ядерной оценки плотности. Оставили. Отдельные ядра. Правильно. Оценка плотности ядра.

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

Определение

Предыдущий рисунок представляет собой графическое представление оценки плотности ядра, которую мы теперь точно определим. Позволять Икс1, Икс2, ..., Иксп быть образец из d-variate случайные векторы взяты из общего распределения, описанного функция плотности ƒ. Оценка плотности ядра определяется как

куда

  • Икс = (Икс1, Икс2, …, Иксd)Т, Икся = (Икся1, Икся2, …, Икся бы)Т, я = 1, 2, …, п находятся d-векторы;
  • ЧАС это полоса пропускания (или сглаживание) d × d матрица, которая симметричный и положительно определенный;
  • K это ядро функция, которая является симметричной многомерной плотностью;
  • .

Выбор функции ядра K не имеет решающего значения для точности оценок плотности ядра, поэтому мы используем стандартную многомерный нормальный ядро во всем: , где H играет роль ковариационная матрица. С другой стороны, выбор матрицы пропускной способности ЧАС является единственным наиболее важным фактором, влияющим на его точность, поскольку он контролирует величину и ориентацию индуцированного сглаживания.[5]:36–39 То, что матрица полосы пропускания также индуцирует ориентацию, является основным отличием многомерной ядерной оценки плотности от ее одномерного аналога, поскольку ориентация не определена для одномерных ядер. Это приводит к выбору параметризации этой матрицы ширины полосы. Три основных класса параметризации (в порядке возрастания сложности): S, класс положительных скаляров, умноженный на единичную матрицу; D, диагональные матрицы с положительными элементами на главной диагонали; и F, симметричные положительно определенные матрицы. В S ядра классов имеют одинаковое сглаживание во всех направлениях координат, D ядра позволяют различное количество сглаживания в каждой из координат, и F ядра позволяют произвольное количество и ориентацию сглаживания. Исторически S и D ядра являются наиболее распространенными из-за вычислительных причин, но исследования показывают, что значительное повышение точности может быть получено с использованием более общего F ядра классов.[6][7]

Comparison of the three main bandwidth matrix parametrisation classes. Оставили. S positive scalar times the identity matrix. Центр. D diagonal matrix with positive entries on the main diagonal. Правильно. F symmetric positive definite matrix.
Сравнение трех основных классов параметризации матрицы пропускной способности. Оставили. S положительный скаляр, умноженный на единичную матрицу. Центр. D диагональная матрица с положительными элементами на главной диагонали. Правильно. F симметричная положительно определенная матрица.

Выбор оптимальной матрицы пропускной способности

Наиболее часто используемый критерий оптимальности для выбора матрицы пропускной способности - это MISE или среднеквадратичная ошибка

Это вообще не имеет выражение в закрытой форме, поэтому обычно используется его асимптотическое приближение (AMISE) в качестве прокси

куда

  • , с р(K) = (4π)−d/2 когда K это нормальное ядро
  • ,
с яd будучи d × d единичная матрица, с м2 = 1 для нормального ядра
  • D2ƒ это d × d Матрица Гессе частных производных второго порядка от ƒ
  • это d2 × d2 матрица интегральных частных производных четвертого порядка от ƒ
  • vec - это векторный оператор, который складывает столбцы матрицы в один вектор, например.

Качество приближения AMISE к MISE[5]:97 дан кем-то

куда о указывает на обычный строчная нотация. Эвристически это утверждение подразумевает, что AMISE является «хорошим» приближением MISE как размера выборки. п → ∞.

Можно показать, что любой разумный селектор полосы пропускания ЧАС имеет ЧАС = О(п−2/(d+4)) где нотация большой O применяется поэлементно. Подставляя это в формулу MISE, получаем, что оптимальным MISE является О(п−4/(d+4)).[5]:99–100 Таким образом, как п → ∞, MISE → 0, т.е. оценка плотности ядра сходится в среднем квадрате а значит, и по вероятности истинной плотности ж. Эти способы сходимости являются подтверждением утверждения в разделе мотивации, что ядерные методы приводят к разумным оценкам плотности. Идеальный селектор оптимальной полосы пропускания - это

Поскольку этот идеальный селектор содержит неизвестную функцию плотности ƒ, его нельзя использовать напрямую. Множество различных разновидностей селекторов полосы пропускания на основе данных возникает из разных оценок AMISE. Мы концентрируемся на двух классах селекторов, которые, как было показано, наиболее широко применимы на практике: сглаженная перекрестная проверка и селекторы плагинов.

Плагин

Подключаемый модуль (PI) оценка AMISE формируется заменой Ψ4 по его оценке

куда . Таким образом это селектор подключаемых модулей.[8][9] Эти ссылки также содержат алгоритмы оптимальной оценки матрицы полосы пропускания пилот-сигнала. грамм и установить, что сходится по вероятности к ЧАСAMISE.

Сглаженная перекрестная проверка

Сглаженная перекрестная проверка (SCV) - это подмножество более крупного класса перекрестная проверка техники. Оценщик SCV отличается от модуля оценщика во втором члене.

Таким образом это селектор SCV.[9][10]Эти ссылки также содержат алгоритмы оптимальной оценки матрицы полосы пропускания пилот-сигнала. грамм и установить, что сходится по вероятности к ЧАСAMISE.

Практическое правило

Эмпирическое правило Сильвермана предлагает использовать куда - стандартное отклонение i-й переменной и . Правило Скотта .

Асимптотический анализ

В разделе выбора оптимальной полосы пропускания мы представили MISE. Его конструкция опирается на ожидаемое значение и отклонение оценщика плотности[5]:97

где свертка оператор между двумя функциями и

Чтобы эти два выражения были четко определены, мы требуем, чтобы все элементы ЧАС стремятся к 0 и что п−1 |ЧАС|−1/2 стремится к 0 как п стремится к бесконечности. Предполагая эти два условия, мы видим, что ожидаемое значение стремится к истинной плотности ж т.е. оценка плотности ядра асимптотически беспристрастный; и что дисперсия стремится к нулю. Использование стандартного разложения среднеквадратичного значения

у нас есть, что MSE стремится к 0, подразумевая, что оценка плотности ядра (среднеквадратичная) согласована и, следовательно, сходится по вероятности к истинной плотности ж. Скорость сходимости MSE к 0 обязательно такая же, как и скорость MISE, отмеченная ранее. О(п−4 / (d + 4)), следовательно, скорость покрытия оценки плотности до ж является Оп(п−2/(d+4)) куда Оп обозначает порядок вероятности. Это устанавливает поточечную сходимость. Функциональное покрытие устанавливается аналогичным образом, рассматривая поведение MISE и отмечая, что при достаточной регулярности интегрирование не влияет на скорость сходимости.

Для рассматриваемых селекторов полосы пропускания на основе данных целью является матрица полосы пропускания AMISE. Мы говорим, что селектор на основе данных сходится к селектору AMISE с относительной скоростью Оп(пα), α > 0, если

Было установлено, что селекторы подключаемого модуля и сглаженной перекрестной проверки (при одной полосе пропускания пилот-сигнала грамм) оба сходятся с относительной скоростью Оп(п−2/(d+6)) [9][11] то есть оба этих основанных на данных селектора являются согласованными оценками.

Оценка плотности с полной матрицей полосы пропускания

Old Faithful Geyser data kernel density estimate with plug-in bandwidth matrix.
Оценка плотности ядра данных Old Faithful Geyser с помощью матрицы пропускной способности плагина.

В пакет ks[12] в р реализует плагин и селекторы сглаженной перекрестной проверки (среди прочего). Этот набор данных (включенный в базовое распределение R) содержит 272 записи с двумя измерениями в каждой: продолжительность извержения (минуты) и время ожидания до следующего извержения (минуты) Старый верный гейзер в национальном парке Йеллоустоун, США.

Фрагмент кода вычисляет оценку плотности ядра с помощью матрицы пропускной способности плагина Опять же, цветные контуры соответствуют наименьшей области, которая содержит соответствующую вероятностную массу: красный = 25%, оранжевый + красный = 50%, желтый + оранжевый + красный = 75%. Чтобы вычислить селектор SCV, Hpi заменяется на Hscv. Это не отображается здесь, поскольку в основном похоже на оценку плагина для этого примера.

библиотека(кс)данные(верный)ЧАС <- Hpi(Икс=верный)фат <- kde(Икс=верный, ЧАС=ЧАС)участок(фат, отображать="fill.contour2")точки(верный, cex=0.5, pch=16)

Оценка плотности с помощью диагональной матрицы ширины полосы

Kernel density estimate with diagonal bandwidth for synthetic normal mixture data.
Оценка плотности ядра с диагональной полосой пропускания для данных синтетической нормальной смеси.

Рассмотрим оценку плотности гауссовой смеси(4π)−1 ехр (-12 (Икс12 + Икс22))+ (4π)−1 ехр (-12 ((Икс1 - 3.5)2 + Икс22)), из 500 случайно сгенерированных точек. Мы используем процедуру Matlab для2-мерные данные.Программа представляет собой метод автоматического выбора полосы пропускания, специально разработанный для гауссовского ядра второго порядка.[13]На рисунке показана оценка плотности стыков, полученная в результате использования автоматически выбранной полосы пропускания.

Скрипт Matlab для примера

Введите следующие команды в Matlab послескачивание и сохранение функции kde2d.min в текущем каталоге.

  Чисто все   % генерировать синтетические данные  данные=[Randn(500,2);      Randn(500,1)+3.5, Randn(500,1);];  % вызвать подпрограмму, которая была сохранена в текущем каталоге   [пропускная способность,плотность,Икс,Y]=kde2d(данные);  % построить данные и оценку плотности  contour3(Икс,Y,плотность,50), держать на  участок(данные(:,1),данные(:,2),'р.','MarkerSize',5)

Альтернативные критерии оптимальности

MISE - это ожидаемый интегрированный L2 расстояние между оценкой плотности и истинной функцией плотности ж. Он является наиболее широко используемым, в основном из-за его управляемости и большей части программного обеспечения, реализующего селекторы полосы пропускания на основе MISE. Существуют альтернативные критерии оптимальности, которые пытаются охватить случаи, когда MISE не является подходящей мерой.[3]:34–37,78 Эквивалент L1 мера, Средняя интегрированная абсолютная ошибка, равна

Его математический анализ значительно сложнее, чем MISE. На практике выигрыш кажется незначительным.[14] В L норма - средняя равномерная абсолютная ошибка

который был исследован лишь кратко.[15] Критерии вероятности ошибки включают критерии, основанные на среднем значении Дивергенция Кульбака – Лейблера

и среднее Расстояние Хеллингера

KL можно оценить с помощью метода перекрестной проверки, хотя селекторы перекрестной проверки KL могут быть неоптимальными, даже если он остается последовательный для ограниченных функций плотности.[16] Селекторы MH были кратко рассмотрены в литературе.[17]

Все эти критерии оптимальности основаны на расстоянии и не всегда соответствуют более интуитивным представлениям о близости, поэтому в ответ на это беспокойство были разработаны более визуальные критерии.[18]

Объективный и управляемый данными выбор ядра

An x-shaped region of empirical characteristic function in Fourier space.
Демонстрация функции фильтра . Квадрат эмпирической функции распределения из N= 10 000 отсчетов «переходного распределения», обсуждаемого в разделе 3.2 (и показанного на рис. 4), для . На этом рисунке представлены две цветовые схемы. Преимущественно темная многоцветная «Х-образная» область в центре соответствует значениям для самого нижнего непрерывного гиперобъема (область, содержащая начало координат); шкала цветов справа относится к цветам в этой области. Слегка окрашенные монотонные области вдали от первого смежного гиперобъема соответствуют дополнительным смежным гиперобъемам (областям) с . Цвета этих областей произвольны и служат только для визуального отличия соседних смежных областей друг от друга.

Недавние исследования показали, что ядро ​​и его пропускная способность могут быть оптимально и объективно выбраны из самих входных данных без каких-либо предположений о форме распределения.[19] Полученная оценка плотности ядра быстро сходится к истинному распределению вероятностей по мере добавления выборок: со скоростью, близкой к ожидается для параметрических оценщиков.[19][20][21] Этот оценщик ядра работает как для одномерных, так и для многомерных выборок. Оптимальное ядро ​​определяется в пространстве Фурье как оптимальная функция демпфирования (преобразование Фурье ядра ) - в терминах преобразования Фурье данных , то эмпирическая характеристическая функция (видеть Оценка плотности ядра ):

[21]

куда, N это количество точек данных, d - количество измерений (переменных), а - фильтр, равный 1 для «принятых частот» и 0 в противном случае. Существуют различные способы определения этой функции фильтра, и простой способ, который работает для одномерных или многомерных выборок, называется «нижний непрерывный гиперобъемный фильтр»; выбирается таким образом, чтобы единственные принятые частоты были непрерывным подмножеством частот, окружающих начало координат, для которых (видеть [21] для обсуждения этой и других функций фильтра).

Обратите внимание, что прямой расчет эмпирическая характеристическая функция (ECF) является медленным, поскольку по существу включает прямое преобразование Фурье выборок данных. Однако было обнаружено, что ECF можно точно аппроксимировать с помощью неоднородное быстрое преобразование Фурье (nuFFT) метод,[20][21] что увеличивает скорость вычислений на несколько порядков (в зависимости от размерности задачи). Комбинация этого объективного метода KDE и приближения ECF на основе nuFFT получила название fastKDE в литературе.[21]

A demonstration of fastKDE relative to a sample PDF. (a) True PDF, (b) a good representation with fastKDE, and (c) a slightly blurry representation.
Нетривиальная смесь нормальных распределений: (а) лежащая в основе PDF, (б) оценка fastKDE для 1 000 000 выборок и (в) оценка fastKDE для 10 000 выборок.

Смотрите также

Рекомендации

  1. ^ Розенблатт, М. (1956). «Замечания о некоторых непараметрических оценках функции плотности». Анналы математической статистики. 27 (3): 832–837. Дои:10.1214 / aoms / 1177728190.
  2. ^ Парзен, Э. (1962). «Об оценке функции плотности вероятности и моды». Анналы математической статистики. 33 (3): 1065–1076. Дои:10.1214 / aoms / 1177704472.
  3. ^ а б Симонов, Дж. (1996). Методы сглаживания в статистике. Springer. ISBN  978-0-387-94716-7.
  4. ^ Сильверман, Б. (1986). Оценка плотности для статистики и анализа данных. Чепмен и Холл / CRC. стр.7–11. ISBN  978-0-412-24620-3.
  5. ^ а б c d Жезл, М.П .; Джонс, М. (1995). Сглаживание ядра. Лондон: Chapman & Hall / CRC. ISBN  978-0-412-55270-0.
  6. ^ Wand, M.P .; Джонс, М. (1993). «Сравнение сглаживающих параметризаций при двумерной оценке плотности ядра». Журнал Американской статистической ассоциации. 88 (422): 520–528. Дои:10.1080/01621459.1993.10476303. JSTOR  2290332.
  7. ^ Duong, T .; Hazelton, M.L. (2003). «Плагин матрицы пропускной способности для двумерной оценки плотности ядра». Журнал непараметрической статистики. 15: 17–30. Дои:10.1080/10485250306039.
  8. ^ Wand, M.P .; Джонс, М. (1994). «Многовариантный выбор пропускной способности плагина». Вычислительная статистика. 9: 97–177.
  9. ^ а б c Duong, T .; Hazelton, M.L. (2005). «Матрицы пропускной способности перекрестной проверки для многомерной оценки плотности ядра». Скандинавский статистический журнал. 32 (3): 485–506. Дои:10.1111 / j.1467-9469.2005.00445.x.
  10. ^ Холл, П .; Marron, J .; Парк, Б. (1992). «Сглаженная перекрестная проверка». Теория вероятностей и смежные области. 92: 1–20. Дои:10.1007 / BF01205233.
  11. ^ Duong, T .; Hazelton, M.L. (2005). «Скорости сходимости для неограниченных матричных селекторов полосы пропускания в многомерной оценке плотности ядра». Журнал многомерного анализа. 93 (2): 417–433. Дои:10.1016 / j.jmva.2004.04.004.
  12. ^ Дуонг, Т. (2007). «ks: оценка плотности ядра и дискриминантный анализ ядра в R». Журнал статистического программного обеспечения. 21 (7). Дои:10.18637 / jss.v021.i07.
  13. ^ Ботев, З.И .; Grotowski, J.F .; Крезе, Д. (2010). «Оценка плотности ядра посредством диффузии». Анналы статистики. 38 (5): 2916–2957. arXiv:1011.2602. Дои:10.1214 / 10-AOS799.
  14. ^ Холл, П .; Жезл, М. (1988). "Сведение к минимуму L1 расстояние в непараметрической оценке плотности ». Журнал многомерного анализа. 26: 59–88. Дои:10.1016 / 0047-259X (88) 90073-5.
  15. ^ Cao, R .; Cuevas, A .; Мантейга, W.G. (1994). «Сравнительное исследование нескольких методов сглаживания при оценке плотности». Вычислительная статистика и анализ данных. 17 (2): 153–176. Дои:10.1016 / 0167-9473 (92) 00066-Z.
  16. ^ Холл, П. (1989). «Об оценке потерь Кульбака-Лейблера и плотности». Анналы статистики. 15 (4): 589–605. Дои:10.1214 / aos / 1176350606.
  17. ^ Ahmad, I.A .; Мугдади, А. (2006). «Взвешенное расстояние Хеллингера как критерий ошибки для выбора полосы пропускания при оценке ядра». Журнал непараметрической статистики. 18 (2): 215–226. Дои:10.1080/10485250600712008.
  18. ^ Marron, J.S .; Цыбаков, А. (1996). «Критерии визуальной ошибки для качественного сглаживания». Журнал Американской статистической ассоциации. 90 (430): 499–507. Дои:10.2307/2291060. JSTOR  2291060.
  19. ^ а б Бернаккья, Альберто; Пиголотти, Симона (01.06.2011). «Самосогласованный метод оценки плотности». Журнал Королевского статистического общества, серия B. 73 (3): 407–422. arXiv:0908.3856. Дои:10.1111 / j.1467-9868.2011.00772.x. ISSN  1467-9868.
  20. ^ а б О’Брайен, Трэвис А .; Коллинз, Уильям Д .; Rauscher, Sara A .; Ринглер, Тодд Д. (01.11.2014). «Снижение вычислительных затрат ECF с помощью nuFFT: быстрый и объективный метод оценки плотности вероятности». Вычислительная статистика и анализ данных. 79: 222–234. Дои:10.1016 / j.csda.2014.06.002.
  21. ^ а б c d е О’Брайен, Трэвис А .; Кашинатх, Картик; Кавано, Николас Р .; Коллинз, Уильям Д .; О’Брайен, Джон П. (2016). «Быстрый и объективный метод многомерной оценки плотности ядра: fastKDE» (PDF). Вычислительная статистика и анализ данных. 101: 148–160. Дои:10.1016 / j.csda.2016.02.014.

внешняя ссылка