Scipy：仅从尾部采样对数正态分布

如何解决Scipy：仅从尾部采样对数正态分布

我正在构建一个模拟，它需要从对数正态分布的尾部随机抽取。选择阈值 τ (tau)，结果条件分布由下式给出：

我需要从该条件分布中随机抽样，其中 F(x) 是具有所选 μ (mu) 和 σ (sigma) 的对数正态分布，τ (tau) 由用户设置。

我现在不雅的解决方案只是从对数正态中采样，丢弃 τ (tau) 以下的任何值，直到我获得所需的样本大小。但我相信这可以改进。

感谢您的帮助！

解决方法

一种干净的方法是定义 rv_continuous 的子类，并实现 _cdf。要绘制变量，您可能还需要定义 _ppf 或 _rvs 方法。

最简单的方法可能是利用 Scipy 提供的截断正态分布。

这给出了以下代码，其中 ν (nu) 作为标准高斯分布的变量，τ (tau) 映射到该分布上的 ν₀。此函数返回一个包含 ranCount 对数正态变量的 Numpy 数组：

import  numpy  as  np
from scipy.stats import truncnorm

def getMySamplesScipy(ranCount,mu,sigma,tau):
    nu0 = (math.log(tau) - mu) / sigma              # position of tau on unit Gaussian
    xs = truncnorm.rvs(nu0,np.inf,size=ranCount)  # truncated unit normal samples
    ys = np.exp(mu + sigma * xs)                    # go back to x space
    return ys

如果由于某种原因这不合适，那么一些常用于高斯变量的技巧，例如 Box-Muller 对截断分布不起作用，但我们总是可以诉诸一般原则：{ {3}} 定理。

因此，我们通过转换统一变量来为我们的变量生成累积概率。我们信任 Scipy，使用它的 erf 误差函数的逆函数从我们的概率返回到 x 空间值。

这给出了类似于以下 Python 代码的内容（没有任何优化尝试）：

import  math
import  random
import  numpy  as  np
import  numpy.random  as  nprd
import  scipy.special as  spfn

# using the "Inverse Method":
def getMySamples(ranCount,tau):
    nu0 = (math.log(tau) - mu) / sigma      # position of tau in standard Gaussian curve
    headCP = (1/2) * (1 + spfn.erf(nu0/math.sqrt(2)))
    tailCP = 1.0 - headCP                   # probability of being in the "tail"

    uvs = np.random.uniform(0.0,1.0,ranCount)  # uniform variates
    cps = (headCP + uvs * tailCP)                # Cumulative ProbabilitieS
    nus = (math.sqrt(2)) * spfn.erfinv(2*cps-1)  # positions in standard Gaussian
    xs  = np.exp(mu + sigma * nus)               # go back to x space
    return xs

替代方案：

我们可以利用与 Inverse Transform Sampling 相关的大量材料。

Zdravko Botev 和 Pierre L'Ecuyer 就该主题发表了一篇相对较新的 (2016) Truncated Gaussian distribution。本文提供了指向公开可用的 review paper 的指针。一些材料非常古老，例如 Luc Devroye 1986 年的书：R source code。

例如，一种可能的基于拒绝的方法：如果τ(tau)映射到标准高斯曲线上的ν₀，则单位高斯分布类似于exp(-ν^{2/2）。如果我们写 ν = ν₀ + δ，这与：exp(-δ²/2) * exp(-ν₀ *δ).}

这个想法是通过参数 ν₀ 的 指数 1 来近似超过 ν₀ 的精确分布。请注意，精确分布始终低于近似分布。然后我们可以随机接受概率为 exp(-δ²/2) 的相对便宜的指数变量。

我们可以选择文献中的等效算法。在 Devroye 书中，第 IX 章第 382 页，有一些伪代码：

重复生成独立的指数随机变量 X 和 Y 直到 X²0²*Y

返回 R 0 + X/ν₀

可以像这样编写 Numpy 再现：

def getMySamplesXpRj(rawRanCount,tau):
    nu0 = (math.log(tau) - mu) / sigma      # position of tau in standard Gaussian
    if (nu0 <= 0):
        print("Error: τ (tau) too small in getMySamplesXpRj")
    rnu0 = 1.0 / nu0

    xs  = nprd.exponential(1.0,rawRanCount)   # exponential "raw" variates
    ys  = nprd.exponential(1.0,rawRanCount)

    allSamples = nu0 + (rnu0 * xs)
    boolArray  = (xs*xs - 2*nu0*nu0*ys) <= 0.0
    samples    = allSamples[boolArray]

    ys  = np.exp(mu + sigma * samples)               # go back to x space
    return ys

根据 Botev-L'Ecuyer 论文中的表 3，该算法的拒绝率非常低。

此外，如果您愿意考虑一些复杂性，还有一些关于用于截断高斯分布的 Non-Uniform Random Variate Generation 的文献，例如 Nicolas Chopin 在 ENSAE-CREST 的 2012 Ziggurat algorithm .

附注：在最新版本的 Python 中，您似乎可以直接使用希腊字母作为变量名，σ 代替 sigma，τ 代替 tau，就像在统计书籍中一样：

$ python3
Python 3.9.6 (default,Jun 29 2021,00:00:00)
>>> 
>>> σ = 2
>>> τ = 7
>>> 
>>> στ = σ * τ
>>> 
>>> στ + 1
15
>>>