postgresql - PostgreSQL 中的 Beta 和 lognorm 分布？

Question

我目前在代码中运行一个相当大的蒙特卡洛模拟，性能还有一些不足之处。

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

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

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

我该怎么做？

Answer 1

方法

正如您所指出的，PostgreSQL 能够生成 Uniform使用 random() 分发功能。

此类问题的一般答案是 Inverse Transform Sampling .
这种方法的局限性是:

显式构造 Quantile Function 的能力(又名 PPF)，可以定义为 Inverse Function的 Improper Integral : PPF(u) = CDF^(-1)(u) | u = CDF(x) = int(PDF(x), x=(-infinty,x)) ;

PostgreSQL mathematical functions的存在|需要构造分位数函数。

也就是说，如果分位数函数是显式的并且我们能够用 PostgreSQL 数学函数构造它，那么我们可以创建一个 Pseudo Random Generator对于使用 random() 的特定分布作为统一 PRG。

简单示例:指数

逆变换采样适用于 Exponential Distribution :

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 的指数分布中抽取的样本.

对数正态

对于 Lognormal distribution分位数函数依赖于 Error Function这在 PostgreSQL 中没有实现。因此我们需要实现缺失的函数(这是一个整体，使用 WINDOWING functions 并非不可能，但可能不是最好的想法)或找到其他方法。

幸运的是，我们可以生成 Normal Distribution使用 Box-Muller Transform 的示例:

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 Distribution因为分位数函数不能表达为一个简单的函数:它需要得到 Regularized Incomplete Beta Function 的倒数。 .我不知道是否有可能:维基百科没有为 Beta 发行版引用分位数函数。

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

技术考虑

正如@Nick Barnes 在他的评论中指出的那样:

使用 random() 的函数不是 IMMUTABLE (它们是 VOLATILE 默认值)因为它们改变了 PostgreSQL PRG 的种子值；

此处介绍的当前实现是幼稚的，它们不处理边缘情况，例如 ln(0.) ;

LANGUAGE SQL 中的函数通常表现良好(尽管我们必须考虑它们的复杂性和收敛性)；

返回 SETOF FLOAT比使用 FLOAT[] 更好并避免需要 unnest() ，就像我在以前版本的 SQL 函数中所做的那样；

限制转换，例如 ::FLOAT只要有可能；

有一个功能pi()无需使用 2.*acos(0.) 对其进行评估.