R - моделировать данные для распределения плотности вероятности, полученные из оценки плотности ядра

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

Итак, я пытаюсь сделать следующее. Я хочу смоделировать данные для всего 250 000 наблюдений, присвоив непрерывное (нецелочисленное) значение в соответствии с оценкой плотности ядра, полученной из эмпирических данных (дискретных), с исходными значениями в диапазоне от -5 до +5. Вот график распределения, который я хочу использовать.

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

Для меня очень важно, чтобы я моделировал новые данные не на основе дискретных вероятностей, а на основе непрерывных, поскольку действительно важно, чтобы значение могло быть, скажем, 2,89, а не 3 или 2. Таким образом, новые значения будут назначаться на основе вероятности, изображенные на сюжете. Наиболее частое значение в смоделированных данных будет где-то около +2, тогда как значения около -4 и +5 будут довольно редкими.

Я довольно много читал о моделировании данных в R и о том, как работают оценки плотности ядра, но на самом деле я вообще не двигаюсь вперед. Итак, мой вопрос в основном состоит из двух шагов - как мне даже смоделировать данные (1) и, кроме того, как мне смоделировать данные, используя это конкретное распределение вероятностей (2)?

Заранее спасибо, я надеюсь, что вы, ребята, можете мне с этим помочь.


person nikUoM    schedule 26.09.2016    source источник
comment
Если я понимаю ваш вопрос, вы, вероятно, могли бы использовать моделирование Монте-Карло. Поскольку ваше распределение известно, вы можете использовать случайную выборку из этого распределения для создания моделируемых данных. Это похоже на то, что вы пытаетесь сделать?   -  person Lloyd Christmas    schedule 26.09.2016
comment
Привет, @LloydChristmas, спасибо за ответ. Дело в том, что распределение известно, но я не хочу моделировать данные, используя фактическое дискретное распределение, а скорее непрерывное, изображенное в ядре. Если это то, о чем вы имеете в виду, это было бы здорово. Как мне это сделать?   -  person nikUoM    schedule 26.09.2016
comment
Здесь есть интересное обсуждение: stats.stackexchange.com/questions/30303/   -  person Lloyd Christmas    schedule 26.09.2016
comment
Я собирался порекомендовать метод отклонения, но он описан в приведенной выше ссылке. Кроме того, @gung предоставил решение, которое может сработать и для вас.   -  person Lloyd Christmas    schedule 26.09.2016
comment
Привет. Я прочитал подход @gung, и я не совсем понимаю его большую часть, но думаю, что это может сработать. Я прочитаю статью, которую он рекомендовал, и посмотрю, смогу ли я создать что-то значимое, используя этот подход. Большое тебе спасибо!   -  person nikUoM    schedule 26.09.2016


Ответы (1)


Используя свои базовые дискретные данные, создайте оценку плотности ядра на столь мелкой сетке, как вы хотите (т. Е. Настолько "близкой к непрерывной", насколько это необходимо для вашего приложения (конечно, в пределах машинной точности и времени вычислений)). Затем выполните выборку из этой плотности ядра, используя значения плотности, чтобы обеспечить более высокую вероятность выборки более вероятных значений вашего распределения. Например:

Поддельные данные, просто чтобы было с чем работать в этом примере:

set.seed(4396)
dat = round(rnorm(1000,100,10))

Создайте оценку плотности ядра. Увеличьте n, если вы хотите, чтобы плотность оценивалась по более мелкой сетке точек:

dens = density(dat, n=2^14)

В этом случае плотность оценивается на сетке из 2 ^ 14 точек с расстоянием mean(diff(dens$x)) = 0,0045 между каждой точкой.

Теперь, выборка из оценки плотности ядра: мы выбираем x-значения оценки плотности и устанавливаем prob равным y-значениям (плотностям) оценки плотности, так что более вероятные x-значения будут с большей вероятностью отобранный:

kern.samp = sample(dens$x, 250000, replace=TRUE, prob=dens$y)

Сравните dens (оценка плотности наших исходных данных) (черная линия) с плотностью kern.samp (красная):

plot(dens, lwd=2)
lines(density(kern.samp), col="red",lwd=2)

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

С помощью описанного выше метода вы можете создать более мелкую сетку для оценки плотности, но вы по-прежнему будете ограничены значениями плотности в точках сетки, используемых для оценки плотности (т. Е. Значениями dens$x). Однако, если вам действительно нужно получить плотность для любого значения данных, вы можете создать функцию аппроксимации. В этом случае вы все равно создадите оценку плотности - при любой ширине полосы и размере сетки, необходимой для захвата структуры данных, - а затем создадите функцию, которая интерполирует плотность между точками сетки. Например:

dens = density(dat, n=2^14)

dens.func = approxfun(dens)

x = c(72.4588, 86.94, 101.1058301)

dens.func(x)
[1] 0.001689885 0.017292405 0.040875436

Вы можете использовать это, чтобы получить распределение плотности при любом значении x (а не только в точках сетки, используемых функцией density), а затем использовать выходные данные dens.func в качестве prob аргумента для sample.

person eipi10    schedule 26.09.2016
comment
Большое спасибо. Я собираюсь попробовать это и посмотреть, что из этого получится. Я вернусь к вам с моими результатами позже. - person nikUoM; 26.09.2016
comment
Уважаемый @ eipi10 - это сработало. Огромное спасибо! - person nikUoM; 27.09.2016