我想从 R 中的自定义分布中提取 1000 个样本
我有以下自定义分布
library(gamlss)
mu <- 1
sigma <- 2
tau <- 3
kappa <- 3
rate <- 1
Rmax <- 20
x <- seq(1, 2e1, 0.01)
points <- Rmax * dexGAUS(x, mu = mu, sigma = sigma, nu = tau) * pgamma(x, shape = kappa, rate = rate)
plot(points ~ x)
如何通过蒙特卡罗模拟从这个分布中随机抽样？
我的第一次尝试是下面的代码，它产生了一个我没想到的直方图形状。
hist(sample(points, 1000), breaks = 51)
这不是我想要的，因为它不遵循与 pdf 相同的分布。
【问题讨论】：
标签： r distribution montecarlo
如果您想要进行蒙特卡罗模拟，您需要从分布中多次抽样，而不是一次抽取大量样本。
您的对象points 的值随着索引增加到400 附近的阈值而增加，然后趋于平稳，然后减小。这就是plot(points ~ x) 显示的内容。它可能描述了一个分布，但points 中值的实际分布是不同的。这显示了值在某个范围内的频率。您会注意到直方图的 x 轴与 plot(points ~ x) 图的 y 轴相似。 points 对象中值的实际分布很容易看到，它与您在随机采样 1000 个值时看到的相似，无需从其中包含 1900 值的对象替换。这是points 中的值分布（无需模拟）：
hist(points, 100)
我故意使用了 100 次休息，以便您可以看到一些细节。
请注意顶部尾部的小凸起，如果您希望直方图看起来像值与索引（或一些增加的 x）的图，您可能不会想到。这意味着points 中有更多的值在2 附近，然后在1 附近。看看你是否可以看看当值在2 附近时plot(points ~ x) 的曲线是如何变平的，以及在0.5 和1.5 之间它是如何非常陡峭的。还要注意直方图低端的大驼峰，并再次查看plot(points ~ x) 曲线。您是否看到大多数值（无论它们是在该曲线的低端还是高端）接近0，或至少小于0.25。如果您查看这些细节，您可能会说服自己直方图实际上正是您所期望的:)
如果您想对该对象的样本进行蒙特卡罗模拟，您可以尝试以下方法：
samples <- replicate(1000, sample(points, 100, replace = TRUE))
如果您想使用points 作为概率密度函数来生成数据，该问题已被提出并回答了here
【讨论】：
points 的直方图是我所期望的。感谢您的出色回答。我现在意识到我要模拟的不是points 的直方图，而是plot(points~x) 的直方图
hist(sample(points, 1000), breaks = 51) 是什么，而不是对实际问题的回答——如何从给定的分布中抽取样本。
让我们将您的（未归一化的）概率密度函数定义为一个函数：
library(gamlss)
fun <- function(x, mu = 1, sigma = 2, tau = 3, kappa = 3, rate = 1, Rmax = 20)
Rmax * dexGAUS(x, mu = mu, sigma = sigma, nu = tau) *
pgamma(x, shape = kappa, rate = rate)
现在一种方法是使用一些MCMC（马尔可夫链蒙特卡罗）方法。例如，
simMCMC <- function(N, init, fun, ...) {
out <- numeric(N)
out[1] <- init
for(i in 2:N) {
pr <- out[i - 1] + rnorm(1, ...)
r <- fun(pr) / fun(out[i - 1])
out[i] <- ifelse(runif(1) < r, pr, out[i - 1])
}
out
}
它从点init 开始，并给N 绘制。该方法可以在许多方面进行改进，但我只是从init = 5 开始，包括 20000 的老化期，并选择每第二次抽奖以减少重复次数：
d <- tail(simMCMC(20000 + 2000, init = 5, fun = fun), 2000)[c(TRUE, FALSE)]
plot(density(d))
【讨论】：
你反转分布的 ECDF：
ecd.points <- ecdf(points)
invecdfpts <- with( environment(ecd.points), approxfun(y,x) )
samp.inv.ecd <- function(n=100) invecdfpts( runif(n) )
plot(density (samp.inv.ecd(100) ) )
plot(density(points) )
png(); layout(matrix(1:2,1)); plot(density (samp.inv.ecd(100) ),main="The Sample" )
plot(density(points) , main="The Original"); dev.off()
【讨论】：
这是另一种方法，它借鉴了R: Generate data from a probability density distribution 和How to create a distribution function in R?：
x <- seq(1, 2e1, 0.01)
points <- 20*dexGAUS(x,mu=1,sigma=2,nu=3)*pgamma(x,shape=3,rate=1)
f <- function (x) (20*dexGAUS(x,mu=1,sigma=2,nu=3)*pgamma(x,shape=3,rate=1))
C <- integrate(f,-Inf,Inf)
> C$value
[1] 11.50361
# normalize by C$value
f <- function (x)
(20*dexGAUS(x,mu=1,sigma=2,nu=3)*pgamma(x,shape=3,rate=1)/11.50361)
random.points <- approx(cumsum(pdf$y)/sum(pdf$y),pdf$x,runif(10000))$y
hist(random.points,1000)
hist((random.points*40),1000) 将像您的原始函数一样获得缩放。
【讨论】：