从R中未知密度的分位数生成随机样本

怎样才能从未知密度的分位数随机抽样数据f(x)对x之间0和4R中？

     f = function(x) ((x-1)^2) * exp(-(x^3/3-2*x^2/2+x))

jlhoward.. 5

如果我正确理解（??），则希望生成随机样本，其分布的密度函数由给出f(x)。一种方法是从均匀分布生成随机样本U[0,1]，然后将该样本转换为您的密度。这是使用cdf的倒数完成的f，该方法之前已在此处进行了描述。

所以让

f(x)     = your density function, 
F(x)     = cdf of f(x), and 
F.inv(y) = inverse cdf of f(x).

在R代码中：

f <- function(x) {((x-1)^2) * exp(-(x^3/3-2*x^2/2+x))}
F <- function(x) {integrate(f,0,x)$value}
F <- Vectorize(F)

F.inv <- function(y){uniroot(function(x){F(x)-y},interval=c(0,10))$root}
F.inv <- Vectorize(F.inv)

x <- seq(0,5,length.out=1000)
y <- seq(0,1,length.out=1000)

par(mfrow=c(1,3))
plot(x,f(x),type="l",main="f(x)")
plot(x,F(x),type="l",main="CDF of f(x)")
plot(y,F.inv(y),type="l",main="Inverse CDF of f(x)")

在上面的代码中，由于f(x)仅在上定义[0,Inf]，我们将其计算F(x)为f(x)0到x 的整数。然后，使用上的uniroot(...)函数将其求反F-y。使用的Vectorize(...)，因为，与几乎所有的R功能，需要integrate(...)与uniroot(...)不上的向量操作。您应该在这些功能上查找帮助文件以获取更多信息。

现在我们只是生成一个随机样本X，U[0,1]并用Z = F.inv(X)

X <- runif(1000,0,1)   # random sample from U[0,1]
Z <- F.inv(X)

最后，我们证明Z确实是作为分布的f(x)。

par(mfrow=c(1,2))
plot(x,f(x),type="l",main="Density function")
hist(Z, breaks=20, xlim=c(0,5))