来自 R 中的自定义分布的模拟答案

【问题标题】：Simulation from a customised distribution in R来自 R 中的自定义分布的模拟
【发布时间】：2017-07-04 11:12:48
【问题描述】：

有谁知道如何从以下概率密度函数（pdf）模拟随机变量：

其中 g(x) 是标准正态分布 N(0,1) 的 pdf。我知道从自定义分布模拟随机变量涉及找到累积密度函数和使用均匀分布（因为在这个平台上也有类似的问题）。但是，我这里的 pdf 看起来比我之前遇到的其他示例要复杂一些。例如，我想找到累积密度函数（需要积分）并不是那么简单。谁能给我建议如何解决这个问题？或者是否有另一种更简单的方法可以从上述分布中进行模拟？期待您的建议。谢谢！

【问题讨论】：

一种方法是设置一个向量 x，其值从（比如说）-5 到 +5 以 0.001 为间隔，计算每个这些值的 f(x)，然后使用 @ 987654325@ 从x 采样，使用f(x) 作为概率。我认为sample 会自动为您标准化概率，以便它们总和为 1。
刚刚尝试绘制此图，它在 +/- 1 处不连续，是吗？

标签： r

【解决方案1】：

您可以尝试马尔可夫链蒙特卡罗方法，例如 Metropolis-Hastings 算法。 https://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm

【讨论】：

【解决方案2】：

这不行吗？

special <- function(x) {
               if (x <= -1) {
                  val <- dnorm(x) / (1+exp(-x-1))
               } else if ((x > -1) & (x < 1)) {
                  val <- dnorm(x) / 2
               } else {
                  val <- dnorm(x) / (1+exp(-x+1))
               }
               return(val)
           }

special(-2)
[1] 0.01452041

special(-0.5)
[1] 0.1760327

【讨论】：

@Palmy 不确定您是想要 pnorm 还是 dnorm。 seankross.com/notes/dpqr看区别

【解决方案3】：

您可以使用sample...

#define some x values
x <- seq(-5, 5, 0.001)
#define the pdf as above
f <- function(x) {
  if (x <= -1) {
    y <- dnorm(x)/(1+exp(-x-1))
  } else if (x >= 1) {
    y <- dnorm(x)/(1+exp(-x+1))
  } else {
    y <- dnorm(x)/2
  }
  return(y)
}
#calculate values of pdf for all x
fx <- f(x)
#create sample using these probabilities
samp <- sample(x, size=10000, replace=TRUE, prob=fx)

hist(samp, breaks=50)

【讨论】：