Почему stat_density (R; ggplot2) и gaussian_kde (Python; scipy) различаются? - PullRequest
Новая
       31

Почему stat_density (R; ggplot2) и gaussian_kde (Python; scipy) различаются?

1 голос
/ 27 марта 2019

Я пытаюсь создать PDF-оценку на основе KDE для ряда распределений, которые могут не распространяться нормально.

Мне нравится, как stat_density ggplot в R, кажется, распознает каждый нарастающий скачок частоты, но не может повторить это с помощью метода Python scipy-stats-gaussian_kde, который кажется слишком гладким.

Я настроил свой код R следующим образом: 

ggplot(test, aes(x=Val, color = as.factor(Class), group=as.factor(Class))) +
             stat_density(geom='line',kernel='gaussian',bw='nrd0' 
                                                            #nrd0='Silverman'
                                                            ,size=1,position='identity')

R example

А мой код на Python: 

kde = stats.gaussian_kde(data.ravel())
kde.set_bandwidth(bw_method='silverman')

python example

Статистика документов показывает здесь , что nrd0 - это метод Сильвермана для настройки bw.

Исходя из кода выше, я использую те же методы kernals (gaussian) и bandwidth (Silverman).

Может кто-нибудь объяснить, почему результаты так разные?

1 Ответ

5 голосов
/ 27 марта 2019

Кажется, существуют разногласия по поводу того, что такое правило Сильвермана.

scipy docs говорят, что правило Сильвермана реализовано как : 

def silverman_factor(self):
    return power(self.neff*(self.d+2.0)/4.0, -1./(self.d+4))

Где d - количество измерений (1, в вашем случае), а neff - эффективный размер выборки (количество точек, при условии отсутствия весов).Таким образом, пропускная способность scipy составляет (n * 3 / 4) ^ (-1 / 5) (умноженное на стандартное отклонение, рассчитанное другим методом).

В отличие от этого, R's stats документы пакета описывают метод Сильвермана как "0,9 разминимальное значение стандартного отклонения и межквартильный интервал, деленные в 1,34 раза на размер выборки на отрицательную одну пятую степени », что также можно проверить в коде R, набрав bw.nrd0 в консоли: 

function (x) 
{
    if (length(x) < 2L) 
        stop("need at least 2 data points")
    hi <- sd(x)
    if (!(lo <- min(hi, IQR(x)/1.34))) 
        (lo <- hi) || (lo <- abs(x[1L])) || (lo <- 1)
    0.9 * lo * length(x)^(-0.2)
}

Википедия , с другой стороны, дает «Правило большого пальца Сильвермана» в качестве одного из многих возможных имен для оценки: 

1.06 * sigma * n ^ (-1 / 5)

Версия википедии эквивалентна версии scipy.

Все три источника цитируют одну и ту же ссылку: Silverman, BW (1986). Оценка плотности для статистики и анализа данных .Лондон: Чепмен и Холл / CRC.п.48. ISBN 978-0-412-24620-3.Википедия и R конкретно ссылаются на страницу 48, тогда как в документах Сципи не упоминается номер страницы.(Я отправил правку в Википедию, чтобы обновить ссылку на ее страницу до стр.45, см. Ниже.)

Приложение

Я нашел PDF изСсылка Сильвермана.

На странице 45 уравнение 3.28 - это то, что используется в статье Википедии: (4 / 3) ^ (1 / 5) * sigma * n ^ (-1 / 5) ~= 1.06 * sigma * n ^ (-1 / 5).Сципи использует тот же метод, переписывая (3 / 4) ^ (-1 / 5) как эквивалент (4 / 3) ^ (1 / 5).Сильверман описывает этот метод следующим образом:

Хотя (3.28) будет хорошо работать, если популяция действительно нормально распределена, она может несколько смягчиться, если популяция является мультимодальной ... поскольку смесь становится более бимодальной,формула (3.28) будет все больше и больше сглаживаться относительно оптимального выбора параметра сглаживания.

Документы scipy ссылаются на эту слабость , заявляя:

Включает автоматическое определение полосы пропускания.Оценка работает лучше всего для унимодального распределения;бимодальные или мультимодальные распределения имеют тенденцию быть чрезмерно сглаженными.

Статья Сильвермана продолжает мотивировать метод, используемый R и Stata.На странице 48 мы получаем уравнение 3.31: 

h = 0.9 * A * n ^ (-1 / 5)
# A defined on previous page, eqn 3.30
A = min(standard deviation, interquartile range / 1.34)

Сильверман описывает этот метод следующим образом:

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

Так что, похоже, Википедия и Сципи используют простую оценкупредложенный Silverman.R и Stata используют более изысканную версию.

...