PostgreSQL 中的 Beta 和 lognorm 分布？

Question

我目前正在代码中运行相当大的蒙特卡洛模拟，并且性能还有待提高。

我想知道是否有一种方法可以直接在数据库上运行它，我认为性能会好得多。我可以生成随机数，但我没有看到统计分布函数。

已经对我有很大帮助的第一步是：

我有一个参数表，其中每一行都是一个 beta 分布及其所有参数。我想用这些分布参数生成随机值并将它们存储在一个单独的表中（蒙特卡罗模拟表，每次模拟运行一行）。

我该怎么办？

Answer 1

方法

正如您所指出的，PostgreSQL 能够使用该函数生成均匀random()分布。

此类问题的一般答案是逆变换采样。这种方法的局限性是：

• 显式构造分位数函数（ 又名 PPF）的能力，该函数可以定义为不当积分的反函数：；PPF(u) = CDF^(-1)(u) | u = CDF(x) = int(PDF(x), x=(-infinty,x))

• 存在构建分位数函数所需的PostgreSQL 数学函数。

也就是说，如果分位数函数是显式的，并且我们能够使用 PostgreSQL 数学函数构造它，那么我们可以使用统一 PRG 为特定分布创建伪随机生成器random()。

简单的例子：指数

逆变换采样非常适合指数分布：

CREATE OR REPLACE FUNCTION expon(N INTEGER, l FLOAT = 1)
RETURNS SETOF FLOAT AS
$BODY$
SELECT
    -(1/l)*ln(1 - random())
FROM
    generate_series(1, N) AS i;
$BODY$
LANGUAGE SQL;

该函数生成N从参数 的指数分布中抽取的样本l。

对数正态

对于对数正态分布，分位数函数依赖于PostgreSQL 中未实现的误差函数。因此，我们需要实现缺失的函数（这是一个整体，使用WINDOWING 函数并非不可能，但可能不是最好的主意）或找到另一种方法。

幸运的是，我们可以使用Box-Muller 变换生成正态分布样本：

CREATE OR REPLACE FUNCTION norm(N INTEGER, mu FLOAT = 0, sigma FLOAT = 1)
RETURNS SETOF FLOAT AS
$BODY$
SELECT
    sigma*sqrt(-2.*ln(random()))*cos(2*pi()*random()) + mu
FROM
    generate_series(1, N) AS i;
$BODY$
LANGUAGE SQL;

以下调用：

SELECT norm(10000);

给出：

而且MLE 的回报(mu=0.021131501222537274, sigma=1.0042820700537662)还不错，我们可能走在好的轨道上。

然后我们可以取这个函数的指数：

CREATE OR REPLACE FUNCTION lognorm(N INTEGER, mu FLOAT = 0, sigma FLOAT = 1)
RETURNS SETOF FLOAT AS
$BODY$
SELECT
    exp(x)
FROM
    norm(N, mu, sigma) AS x;
$BODY$
LANGUAGE SQL;

我们有一个用于对数正态分布的 PRG。

以下调用：

SELECT lognorm(10000);

也给出了可接受的结果：

MLE 返回(sigma=0.9996878296400589, loc=0.0, exp(mu)=1.0002728392916154)。

数值积分和误差函数

虽然它可能性能不高，但使用Trapezoid Rule来估计 PostgreSQL 的误差函数是相当容易的。认为这是一个幼稚的实现：

CREATE OR REPLACE FUNCTION erf(x FLOAT, dx NUMERIC = 1e-3)
RETURNS FLOAT AS
$BODY$
WITH
D AS (
SELECT
    y::FLOAT,
    exp(-((y::FLOAT)^2)) AS fx0,
    LEAD(exp(-((y::FLOAT)^2))) OVER(ORDER BY y) AS fx1
FROM
    generate_series(0, x::NUMERIC, dx) AS y
)
SELECT
    COALESCE((2/sqrt(pi()))*SUM((D.fx1 + D.fx0)*dx::FLOAT/2), 0.)
FROM D;
$BODY$
LANGUAGE SQL IMMUTABLE;

如果我们将结果与精确形式（Python、scipy）进行比较，看起来还不错，我们至少得到了 6 位有效数字：

      x      psql     scipy        errabs        errrel
0   0.0  0.000000  0.000000  0.000000e+00           NaN
5   0.5  0.520500  0.520500 -7.323189e-08 -1.406953e-07
10  1.0  0.842701  0.842701 -6.918458e-08 -8.209863e-08
15  1.5  0.966105  0.966105 -2.973257e-08 -3.077571e-08
20  2.0  0.995322  0.995322 -6.888995e-09 -6.921371e-09
25  2.5  0.999593  0.999593 -9.076190e-10 -9.079885e-10
30  3.0  0.999978  0.999978 -6.962642e-11 -6.962795e-11
35  3.5  0.999999  0.999999 -3.149592e-12 -3.149594e-12
40  4.0  1.000000  1.000000 -8.404388e-14 -8.404388e-14
45  4.5  1.000000  1.000000  1.110223e-16  1.110223e-16
50  5.0  1.000000  1.000000  2.442491e-15  2.442491e-15

因此，我们可以使用erf函数对正态和对数正态执行逆变换采样，就像我们对指数所做的那样，但我可能是一个坏主意。由于算法复杂性和集成不准确，它的性能应该很差。

贝塔

不幸的是，逆变换采样似乎不适合Beta 分布，因为分位数函数不能表示为简单函数：它需要获得正则化不完全 Beta 函数的逆函数。我不知道这是否可能：维基百科没有为 Beta 分布引用的分位数函数。

对于这种情况，您可能需要用某种编程语言（例如 C/C++）编译该函数，并将其绑定到 PostgreSQL 函数，如 @Nick Barnes 在他的评论中建议的那样。

技术考虑

正如@Nick Barnes 在他的评论中指出的：

• 函数使用random()不是IMMUTABLE（它们是VOLATILE默认的），因为它们改变了 PostgreSQL PRG 的种子值；
• 这里介绍的当前实现很幼稚，它们不处理边缘情况，例如ln(0.)；
• 中的函数LANGUAGE SQL通常表现良好（尽管我们必须考虑它们的复杂性和收敛性）；
• ReturningSETOF FLOAT比 using 更好FLOAT[]，并且避免了需要unnest()，就像我在 SQL 函数的早期版本中所做的那样；
• ::FLOAT尽可能限制演员；
• 有一个函数pi()不需要用 来评估它2.*acos(0.)。