模拟二元帕累托分布

Question

我正在寻找一个包/代码，当两个随机变量相关时（并且相关性可以由用户指定），它会生成二元帕累托分布。非常感谢您的帮助！

Answer 1

我找不到任何现有的 R 包，但发现您的问题很有趣；所以我想向您展示一种如何使用inverse transform sampling 从bivariate type I Pareto distribution 中采样的方法。

---

理论

I 型二元帕累托分布的联合 pdf 由下式给出

这里的目标是

1. 从边际分布f(x2)中为x2抽取样本，然后
2. 从条件分布 f(x1|x2) 给定 x2 为 x1 抽取样本。

边际分布和条件分布由下式给出（参见例如[Mardia, Annals of Mathematical Statistics 33, 1008 (1962)]）

我们可以使用逆变换采样来抽取样本，这需要边缘分布和条件分布的累积分布函数。这很容易计算，我们得到

然后 x1 和 x2 的样本由

给出

其中 u 是区间 [0,1] 内标准均匀分布的随机数。

---

R 实现

1. 我们定义了两个函数来使用上述逆变换采样从边缘分布和条件分布中对 x1 和 x2 的值进行采样。

   rpareto_inv <- function(n, theta, a) {
       u <- runif(n, min = 0, max = 1);
       return(theta / (u ^ (1 / a)));
   }
   
   
   rpareto_cond_inv <- function(x2, theta1, theta2, a) {
     u <- runif(length(x2), min = 0, max = 1);
     return(theta1 + theta1 / theta2 * x2 * (1 / (u ^ (1 / (a + 1))) - 1));
   }
   

2. 我们为采样和分布参数选择了一些值：

   n <- 10^5;     # Number of samples
   theta1 <- 5;   # Location parameter 1
   theta2 <- 2;   # Location parameter 2
   a <- 3;        # Shape parameter
   

3. 现在我们可以绘制样本了

   set.seed(2017);
   x2 <- rpareto_inv(n, theta = theta2, a = a);
   x1 <- rpareto_cond_inv(x2, theta1, theta2, a);
   

4. 我们可以显示二维密度图，并将一些样本汇总统计数据与其理论（总体）值进行比较。

   require(ggplot2);
   df <- cbind.data.frame(x1 = x1, x2 = x2);
   ggplot(df, aes(x1, x2)) +
       geom_density_2d() +
       xlim(theta1, 1.5 * theta1) +
       ylim(theta2, 1.5 * theta2);
   

   metrics <- cbind.data.frame(
       obsrv = c(mean(df$x1), mean(df$x2), cor(df$x1, df$x2), cov(df$x1, df$x2)),
       theor = c(a * theta1 / (a - 1), a * theta2 / (a - 1), 1/a, theta1 * theta2 / ((a - 1)^2 * (a - 2))));
   rownames(metrics) <- c("Mean(x1)", "Mean(x2)", "Correlation", "Covariance")
   #                obsrv     theor
   #Mean(x1)    7.4947124 7.5000000
   #Mean(x2)    3.0029318 3.0000000
   #Correlation 0.3429634 0.3333333
   #Covariance  2.3376545 2.5000000
   

   你可以看到协议很好。另请注意，x1 和 x2 之间的相关性由尺度参数 a 表征。因此，如果您想模拟具有特定相关性 r 的二元帕累托分布的数据，只需将 shape 参数设置为 1/r。有关分布和其他汇总统计的更多详细信息，请参阅[Mardia, Annals of Mathematical Statistics 33, 1008 (1962)]。

最后，您还可以使用简单的接受-拒绝采样方法，但我想这比我在这里展示的逆变换采样方法慢很多。