[论文解读] An efficient algorithm for sampling from $\sin^k(x)$ for generating random correlation matrices

Enes Makalic, Daniel F. Schmidt|arXiv (Cornell University)|Sep 14, 2018

Cellular Automata and Applications参考文献 1被引用 4

一句话总结

本文提出了一种高效的拒绝采样算法，用于在 $(0, \pi)$ 上从与 $\sin^k(x)$ 成比例的分布中生成随机变量，采用贝塔分布作为包络。该方法在所有 $k \geq 1$ 下实现了最大平均接受率 $\pi/(2\sqrt{2}) \approx 1.11$ 次迭代/样本，从而通过超球面参数化实现了高维随机相关矩阵的快速生成。

ABSTRACT

In this note, we develop a novel algorithm for generating random numbers from a distribution with a probability density function proportional to $\sin^k(x)$, $x \in (0,π)$ and $k \geq 1$. Our algorithm is highly efficient and is based on rejection sampling where the envelope distribution is an appropriately chosen beta distribution. An example application illustrating how the new algorithm can be used to generate random correlation matrices is discussed.

研究动机与目标

• 开发一种计算高效的采样方法，用于从 $(0, \pi)$ 上与 $\sin^k(x)$ 成比例的概率密度函数中抽样，该分布出现在随机相关矩阵生成中。
• 克服逆变换采样方法的低效性，后者需要 $O(k)$ 次操作，对于大 $k$ 值变得不切实际。
• 提供一种数值稳定的替代方案，以避免数值积分在大 $k$ 值下出现的精度问题。
• 通过高效采样所需的角变量 $\theta_{ij}$，实现高维相关矩阵的可扩展生成。
• 在 R 和 MATLAB 等广泛使用的平台中实现并基准测试该算法，以支持实际应用。

提出的方法

• 该算法使用拒绝采样，以对称的 $\text{Beta}(k+1, k+1)$ 分布作为包络密度。
• 接受准则基于比值 $f_X(x)/g_Y(x)$，若满足 $\frac{1}{k}\log U \leq \log\left(\frac{\pi^2 \sin X}{4X(\pi - X)}\right)$ 则接受。
• 通过解析推导得出上界 $M_k = \frac{f_X(\pi/2)}{g_Y(\pi/2)}$，确保包络条件 $f_X(x) \leq M_k g_Y(x)$ 成立。
• 使用伽马函数和贝塔函数表达 $f_X(x)$ 和 $g_Y(x)$ 的归一化常数，实现精确计算。
• 通过将该方法应用于 Pourahmadi 和 Wang 提出的超球面参数化方法，对 $p \times p$ 相关矩阵进行生成，从而验证了该方法。
• 该算法已在 R（通过 randcorr 包）和 MATLAB（文件 ID: 68810）中实现，供公众使用。

实验结果

研究问题

• RQ1 能否设计一种拒绝采样方法，以高效地从 $(0, \pi)$ 上的 $\sin^k(x)$ 分布中抽样，适用于所有 $k \geq 1$？
• RQ2 何种包络分布能提供对接受率的紧致上界，同时允许快速抽样？
• RQ3 在所有 $k \geq 1$ 的情况下，理论上最大平均每次采样的迭代次数是多少？
• RQ4 与逆变换采样相比，该方法在高维相关矩阵生成中的性能如何？
• RQ5 该算法能否高效扩展，用于生成 $p$ 最大达 1000 的大 $p \times p$ 相关矩阵？

主要发现

• 最大平均每次采样的迭代次数为 $\pi/(2\sqrt{2}) \approx 1.11$，在 $k \to \infty$ 的极限下达到。
• 当 $k=1$ 时，平均迭代次数为 $\pi/3 \approx 1.047$，并随 $k$ 单调递增。
• 该算法在普通笔记本电脑上生成一个 $100 \times 100$ 的随机相关矩阵耗时约 0.02 秒。
• 生成一个 $1000 \times 1000$ 的相关矩阵耗时约 0.5 秒，表明其与矩阵大小呈线性缩放关系。
• 与逆变换采样相比，该方法性能更优：在相同硬件条件下，后者生成 $100 \times 100$ 矩阵耗时 0.9 秒，生成 $600 \times 600$ 矩阵耗时 42 秒。
• 该拒绝采样器具有数值稳定性，避免了大 $k$ 值下数值积分的精度问题。