对于相关分布抽样，是否有 scipy _norm_pdf 的快速替代方案？

Question

我已经为蒙特卡洛模拟拟合了一系列 SciPy 连续分布，并希望从这些分布中获取大量样本。但是，我希望能够获取相关样本，以便第 i 个样本从每个分布中获取例如第 90 个百分位数。

在这样做的过程中，我发现了 SciPy 性能中的一个怪癖：

# very fast way to many uncorrelated samples of length n
for shape, loc, scale, in distro_props:
    sp.stats.norm.rvs(*shape, loc=loc, scale=scale, size=n)

# verrrrryyyyy slow way to take correlated samples of length n
correlate = np.random.uniform(size=n)
for shape, loc, scale, in distro_props:
    sp.stats.norm.ppf(correlate, *shape, loc=loc, scale=scale)

大多数关于此的结果声称这些 SciPy 发行版的缓慢来自类型检查等包装器。然而，当我分析代码时，大部分时间花在了底层数学函数 [_continuous_distns.py:179(_norm_pdf)]1 上。此外，它按 n 缩放，这意味着它在内部循环遍历每个元素。

SciPy docs on rv_continuous 几乎似乎暗示子类应该覆盖它以提高性能，但我将 monkeypatch 到 SciPy 以加速他们的 ppf 似乎很奇怪。我只想从 ppf 公式计算法线，但我也使用对数正态和偏斜法线，这更难实现。

那么，Python 中计算正态分布、对数正态分布和偏正态分布的快速 ppf 的最佳方法是什么？或者更广泛地说，从几个这样的分布中获取相关样本？

Answer 1

如果你只需要普通的ppf，那么慢确实令人费解，但你可以使用scipy.special.erfinv代替：

x = np.random.uniform(0,1,100)
np.allclose(special.erfinv(2*x-1)*np.sqrt(2),stats.norm().ppf(x))
# True
timeit(lambda:stats.norm().ppf(x),number=1000)
# 0.7717257660115138
timeit(lambda:special.erfinv(2*x-1)*np.sqrt(2),number=1000)
# 0.015020604943856597

编辑：

lognormal 和 triangle 也很简单：

c = np.random.uniform()

np.allclose(np.exp(c*special.erfinv(2*x-1)*np.sqrt(2)),stats.lognorm(c).ppf(x))
# True

np.allclose(((1-np.sqrt(1-(x-c)/((x>c)-c)))*((x>c)-c))+c,stats.triang(c).ppf(x))
# True

偏正可惜我还不够熟

Answer 2

最终，这个问题是我使用 skew-normal 发行版引起的。偏态法线的 ppf 实际上没有封闭形式的解析定义，因此为了计算 ppf，它退回到 scipy.continuous_rv 的数值近似，其中涉及迭代计算 cdf 并使用它来ppf 值归零。 skew-normal pdf是normal pdf和normal cdf的乘积，所以这个数值近似多次调用了normal的pdf和cdf。这就是为什么当我分析代码时，它 看起来 正态分布是问题所在，而不是 SKU 正态分布。这个问题的另一个答案是能够通过跳过类型检查来节省时间，但实际上并没有对运行时间增长产生影响，只是在小 n 运行时间上有所不同.

为了解决这个问题，我用 Johnson SU 分布替换了偏态正态分布。它比正态分布多了 2 个自由参数，因此可以有效地适应不同类型的偏斜和峰态。它支持所有实数，并且具有封闭形式的 ppf 定义，可在 SciPy 中快速实现。下面您可以看到我从第 10、50 和 90 个百分位数拟合的 Johnson SU 分布示例。

对于相关分布抽样，是否有 scipy _norm_pdf 的快速替代方案？

Is there a fast alternative to scipy _norm_pdf for correlated distribution sampling?

python

numpy

scipy

distribution

montecarlo