如何从R中已知的双变量PDF生成随机变量

Question

我在DX x DY矩形区域中有一个二元概率密度函数：

链接到我的pdf

我正在使用R。我如何在这个pdf分布之后的矩形内生成随机(x，y)点？

我读过很多关于"逆变换采样"的答案，但我没有单变量pdf。我见过这种方法，但看起来很乏味，而且实现起来非常困难。

提前谢谢！

Answer 1

这似乎更像是一个统计问题，而不是一个编程问题。无论如何，拒绝抽样应该在你的情况下工作。你可以在维基百科上阅读详细信息。以下函数执行拒绝采样(n…样本数量；pdf…概率密度函数；maxval…pdf可以产生的最大(或更大)值；xlim，ylim。。。密度函数产生大于零的值的边界框)：

reject.sample.2d <- function(n,pdf,maxval,xlim,ylim)
{
smpl <- data.frame(x=numeric(n),y=numeric(n))
i <- 0
while (i<n)
{
xval <- runif(1,xlim[1],xlim[2])
yval <- runif(1,ylim[1],ylim[2])
if (runif(1)<pdf(xval,yval)/maxval)
{
i <- i+1
smpl[i,] <- c(xval,yval)
}
}
return(smpl)
}

例如，它可以用于二维正态分布

mydens <- function(x,y)
{
dnorm(x)*dnorm(y)
}

这样(这里定义的2d正态分布的最大值在x=0，y=0并且比0.16小一点)：

res <- reject.sample.2d(5000,mydens,0.16,c(-5,5),c(-5,5))

一些检查：

> sd(res[["x"]])
[1] 1.015413
> sd(res[["y"]])
[1] 0.9981738
> shapiro.test(res[["x"]])
Shapiro-Wilk normality test
data:  res[["x"]]
W = 0.9995, p-value = 0.1603
> shapiro.test(res[["y"]])
Shapiro-Wilk normality test
data:  res[["y"]]
W = 0.9997, p-value = 0.8304

p.s。也可以使用逆变换采样。计算沿x轴的边际分布。使用此分布进行逆变换采样以获得x值。对于该x值，使用条件概率分布(给定所获得的x的y的概率)，并再次使用逆变换采样，现在用于y值。