加快从叠加截断正态分布中绘制随机值

Posted 2023-03-12

【中文标题】加快从叠加截断正态分布中绘制随机值

【英文标题】：Speed up drawing random values from superimposed truncated normal distributions

【发布时间】：2021-05-05 14:20:24

【问题描述】：

我想从一个分布中抽取 N 个随机样本，该分布是两个截断的正态分布之和。我通过从 scipy.stats 子类化 rv_continuous 类并提供一个 pdf 来获得我想要的东西，该 pdf 是两个给定 pdf 的平均值：

import numpy as np
from scipy import stats

my_lim = [0.05, 7]  # lower and upper limit
my_loc = [1.2, 3]  # loc values of the two truncated normal distributions
my_scale = [0.6, 2]  # scale values of the two truncated normal distributions

class sum_truncnorm(stats.rv_continuous):
    def _pdf(self, x):
        return (stats.truncnorm.pdf(x,
                                    a=(my_lim[0] - my_loc[0]) / my_scale[0],
                                    b=(my_lim[1] - my_loc[0]) / my_scale[0],
                                    loc=my_loc[0], 
                                    scale=my_scale[0]) + 
                stats.truncnorm.pdf(x,
                                    a=(my_lim[0] - my_loc[1]) / my_scale[1],
                                    b=(my_lim[1] - my_loc[1]) / my_scale[1],
                                    loc=my_loc[1],  
                                    scale=my_scale[1]) / 2

但是，使用：

my_dist = sum_truncnorm()
my_rvs = my_dist.rvs(size=10)

非常慢，每个随机值大约需要 5 秒。

我确信这可以更快地完成，但我不知道该怎么做。我是否应该将我的分布定义为（非截断）正态分布的总和并在之后强制截断？我在这个方向上做了一些测试，但这只是快了 10 倍左右，因此仍然很慢。

Google 告诉我，我可能需要使用逆变换采样并覆盖 _rvs 方法，但我未能使其适用于我的截断分布。

【参考方案1】：

首先，您必须确保 _pdf 已标准化。 框架不会检查它，否则会默默地产生废话。

其次，要使绘图快速变化，您需要实现 _ppf 或 _rvs。仅使用 _pdf，它会通过通用代码路径（数字集成和根查找），这就是您当前版本缓慢的原因。

【讨论】：

两个 truncnorm 的 _pdf 应该归一化吧！？当我对 _pdf 使用相同的范围并将两者的总和除以 2 时，我假设它已正确归一化。那不正确？关于_rvs的实现：这就是我尝试但未能实现的。

我怀疑使用docs.scipy.org/doc/scipy/reference/generated/…写下_ppf函数更容易

您使用erf 来获得截断高斯的归一化。