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

Если у меня есть функция плотности и я рисую ее с определенной полосой пропускания, я визуально определяю, что есть 7 локальных максимумов. Я просто хотел бы знать, как строить отдельные распределения конкретных максимумов на одном и том же графике. Кроме того, можно ли точно узнать, где возникают максимумы, запустив какой-либо код? Я могу делать приблизительные оценки, используя график, но есть ли функция R, которую я могу использовать для получения точных точек? Я хотел бы знать среднее значение и дисперсию 7 плотностей, которые я идентифицировал. В частности, у меня есть следующее:

plot(density(stamp, bw=0.0013,kernel = "gaussian"))

r plot statistics kernel-density
person A.Chandy    schedule 23.10.2016    source источник

Ответы (2)


arrow_upward
3
arrow_downward

Определение того, какие режимы являются реальными в оценке плотности ядра, зависит от того, какую полосу пропускания вы выбрали. Это сложная вещь, и я не советую выбирать только одну полосу пропускания, так как даже разные оптимальные эмпирические правила могут дать вам разные ответы. В общем, количество мод kde меньше, чем количество базовой плотности в случае пересглаженного и больше в случае недосглаженного. Есть много статей, посвященных этой теме и предлагающих множество вариантов, помогающих определить достоверность режима. например, ознакомьтесь с тестом режимов Сильвермана для ядер Гаусса, алгоритмом простых чисел Фридмана и Фишера, sZer Маррона и деревом режимов Миннотта и Скотта.

Наивная вещь, которую вы можете сделать, учитывая единственный выбор пропускной способности KDE, - это проверить длину прогона.

На самом деле с выбранной вами полосой пропускания я нахожу 9 режимов. Просто рассчитайте изменение знака разницы в серии и вычислите совокупную длину прогонов, чтобы найти точки. Каждая вторая точка будет модусом или антимодусом, в зависимости от того, что наступило раньше. (Вы можете проверить знак, чтобы определить это)

library(BSDA)
dstamp <- density(Stamp$thickness, bw=0.0013, kernel = "gaussian")
chng <- cumsum(rle(sign(diff(dstamp$y)))$lengths)
plot(dstamp)
abline(v = dstamp$x[chng[seq(1,length(chng),2)]])

введите описание изображения здесь

person shayaa    schedule 23.10.2016
comment
Спасибо за ваш ответ. Как мне узнать средства 9 режимов, которые вы идентифицировали? - person A.Chandy; 23.10.2016
comment
Кроме того, есть ли способ узнать ровно 7 режимов? Я думал, что моя пропускная способность дала мне 7, но оказалось, что это 9 - person A.Chandy; 23.10.2016
comment
вы можете использовать mean и sort соответственно. - person shayaa; 24.10.2016

arrow_upward
0
arrow_downward

Поскольку мне нужно было что-то, чтобы получить только самые сильные моды, я создал очень простой алгоритм, который позволяет вам увеличить чувствительность, настраивая количество выборок плотности (чтобы уменьшить локальный шум) и установить порог минимальной плотности, пропорциональный максимальной плотности (чтобы уменьшить общий шум).

find_posterior_modes <- function(x, n.samples = 100, filter = .1) {
    d <- density(x, n = n.samples)
    x <- with(d, sapply(2:(n.samples - 1), function(i) if (y[i] > y[i - 1] & y[i] > y[i + 1] & y[i] > max(y) * filter) x[i]))
    unlist(x)
}
person Bakaburg    schedule 02.10.2020