Как лучше всего моделировать произвольную одномерную случайную переменную с использованием ее вероятностной функции?

Ответ 1

Вот (медленная) реализация обратного метода cdf, когда вам дается только плотность.

den<-dnorm #replace with your own density

#calculates the cdf by numerical integration
cdf<-function(x) integrate(den,-Inf,x)[[1]]

#inverts the cdf
inverse.cdf<-function(x,cdf,starting.value=0){
 lower.found<-FALSE
 lower<-starting.value
 while(!lower.found){
  if(cdf(lower)>=(x-.000001))
   lower<-lower-(lower-starting.value)^2-1
  else
   lower.found<-TRUE
 }
 upper.found<-FALSE
 upper<-starting.value
 while(!upper.found){
  if(cdf(upper)<=(x+.000001))
   upper<-upper+(upper-starting.value)^2+1
  else
   upper.found<-TRUE
 }
 uniroot(function(y) cdf(y)-x,c(lower,upper))$root
}

#generates 1000 random variables of distribution 'den'
vars<-apply(matrix(runif(1000)),1,function(x) inverse.cdf(x,cdf))
hist(vars)

Ответ 2

Чтобы прояснить "использование Metropolis-Hastings" выше:

предположим, что ddist() - ваша функция плотности вероятности

что-то вроде:

n <- 10000
cand.sd <- 0.1
init <- 0
vals <- numeric(n)
vals[1] <- init 
oldprob <- 0
for (i in 2:n) {
    newval <- rnorm(1,mean=vals[i-1],sd=cand.sd)
    newprob <- ddist(newval)
    if (runif(1)<newprob/oldprob) {
        vals[i] <- newval
    } else vals[i] <- vals[i-1]
   oldprob <- newprob
}

Примечания:

• полностью непроверенный
• эффективность зависит от распределения кандидата (т.е. значения cand.sd). Для максимальной эффективности настройте cand.sd на скорость приема 25-40%
• результаты будут автокоррелированы... (хотя, я думаю, вы всегда могли sample() результаты для их скремблирования или тонкие)
• может потребоваться отменить "ожог", если ваше начальное значение является странным.

Классический подход к этой проблеме - выборка отбраковки (см., например, Press et al. Numerical Recipes)

Ответ 3

Использовать функцию кумулятивного распределения http://en.wikipedia.org/wiki/Cumulative_distribution_function

Тогда просто используйте его обратный. Проверьте здесь лучшее изображение http://en.wikipedia.org/wiki/Normal_distribution

Это означает: выберите случайное число из [0,1] и установите как CDF, затем проверьте значение

Он также называется квантильной функцией.

Ответ 4

Это комментарий, но у меня недостаточно репутации, чтобы оставить комментарий к ответу Бена Болкера.

Я новичок в Метрополис, но ИМХО этот код неверен, потому что:

а) newval взят из нормального распределения, тогда как в других кодах он взят из равномерного распределения; это значение должно быть взято из диапазона, охватываемого случайным числом. Например, для гауссовского распределения это должно быть что-то вроде runif (1, -5, +5).

б) значение вероятности должно обновляться только в случае принятия.

Надеюсь, что это поможет, и надеюсь, что кто-то с репутацией может исправить этот ответ (особенно мой, если я ошибаюсь).

# the distribution 
ddist <- dnorm
# number of random number
n <- 100000
# the center of the range is taken as init
init <- 0
# the following should go into a function
vals <- numeric(n)
vals[1] <- init 
oldprob <- 0
for (i in 2:n) {
  newval <- runif(1, -5, +5)
  newprob <- ddist(newval)
  if (runif(1) < newprob/oldprob) {
    vals[i] <- newval
    oldprob <- newprob
  } else vals[i] <- vals[i-1]
}
# Final view
hist(vals, breaks = 100)
# and comparison
hist(rnorm(length(vals)), breaks = 100)

Ответ 5

Вы можете использовать metropolis-hastings для получения образцов из плотности.