Генерация случайных величин из функции распределения с использованием обратной выборки

У меня есть конкретная функция плотности, и я хочу сгенерировать случайные величины, зная выражение функции плотности.

Например, функция плотности:

df=function(x) { - ((-a1/a2)*exp((x-a3)/a2))/(1+exp((x-a3)/a2))^2 }

Из этого выражения я хочу сгенерировать 1000 случайных элементов с таким же распределением.

Я знаю, что мне следует использовать метод обратной выборки. Для этого я использую функцию CDF моего PDF-файла, которая рассчитывается следующим образом:

cdf=function(x) { 1 - a1/(1+exp((x-a3)/a2))

Идея состоит в том, чтобы сгенерировать равномерно распределенные образцы, а затем сопоставить их с моими функциями CDF, чтобы получить обратное сопоставление. Что-то вроде этого:

random.generator<-function(n) sapply(runif(n),cdf)

а затем вызовите его с желаемым количеством случайных величин для генерации.

random.generator(1000)

Это правильный подход?


r random probability-density
person Lydie    schedule 02.02.2016    source источник
comment
Вам необходимо рассчитать или оценить обратную функцию CDF.   -  person RHertel    schedule 02.02.2016

Ответы (1)


arrow_upward
4
arrow_downward

Первый шаг - сделать обратную функцию вашей cdf, что в этом случае можно сделать с помощью простой арифметики:

invcdf <- function(y) a2 * log(a1/(1-y) - 1) + a3

Теперь вы хотите вызвать обратный cdf со стандартными равномерно распределенными случайными величинами для выборки:

set.seed(144)
a1 <- 1 ; a2 <- 2 ; a3 <- 3
invcdf(runif(10))
#  [1] -2.913663  4.761196  4.955712  3.007925  1.472119  4.138772 -3.568288
#  [8]  4.973643 -1.949684  6.061130

Это гистограмма 10000 смоделированных значений:

hist(invcdf(runif(10000)))

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

А вот и сюжет pdf:

x <- seq(-20, 20, by=.01)
plot(x, df(x))

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

person josliber♦    schedule 02.02.2016
comment
спасибо за подробный ответ, это очень полезно. Я вижу, что вы внесли изменения как в мою функцию cdf, так и в предложенную вами обратную! Есть ли ошибка в моем выражении CDF? - person Lydie; 02.02.2016
comment
чтобы объяснить больше, у меня изначально есть функция CDF, определенная следующим образом: cdf = function (x) {1 - a1 / (1 + exp ((x-a3) / a2)) Функция плотности задается производным: df = function (x) {((-a1 / a2) * exp ((x-a3) / a2)) / (1 + exp ((x-a3) / a2)) ^ 2} - person Lydie; 02.02.2016
comment
@ N.Fk Да, изначально ваш cdf принимал значения, близкие к 1 для низких входов и значения, близкие к 0 для высоких входов, поэтому я думаю, что вы случайно его перевернули - я просто сделал одно за вычетом выражения, которое у вас было раньше. - person josliber♦; 02.02.2016
comment
Если я хочу сделать какие-то условия в сгенерированных числах, знаете ли вы, как это сделать? Я имею в виду, что, например, я хочу, чтобы сгенерированные значения были больше или равны 100? Это возможно ? Я уже пытаюсь повторять процесс до тех пор, пока не получу значения, удовлетворяющие моим ограничениям, но это неэффективно! - person Lydie; 04.02.2016
comment
@ Is.Fk Поскольку это звучит как другой вопрос, могу ли я предложить вам задать его отдельно, используя кнопку «Задать вопрос»? Таким образом, все сообщество сможет ответить на него. - person josliber♦; 04.02.2016