Scipy：仅来自尾部的样本对数正态分布

Question

我正在建立一个模拟，需要从对数正态分布的尾部随机抽取。选择阈值τ（tau），得到的条件分布如下：

我需要从条件分布中随机取样，其中F（x）是对数正态分布，选择µ（mu）和σ（sigma），τ（tau）由用户设置

我现在的不雅观的解决方案就是从对数正态分布中取样，去掉τ（tau）下的任何值，直到得到我需要的样本量。但我相信这是可以改进的

谢谢你的帮助

Answer 1

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

Answer 2

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

这给出了以下代码，其中ν（nu）作为标准高斯分布的变量，τ（tau）映射到该分布上的ν₀。此函数返回包含ranCount lognormal变量的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对截断分布不起作用，但是我们可以始终求助于一个一般原则：Inverse Transform Sampling定理

所以我们通过变换均匀变量来生成变量的累积概率。我们信任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, mu, sigma, 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

备选方案：

我们可以利用大量与Truncated Gaussian distribution相关的材料

Zdravko Botev和Pierre L'Ecuyer最近（2016年）就这一主题发表了一篇文章。本文提供了一个指向公共可用的R source code的指针。有些材料非常陈旧，例如1986年吕克·德夫罗耶的书：Non-Uniform Random Variate Generation

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

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

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

重复 生成独立的指数随机变量X和Y 直到X²<；=2*ν₀²*Y

返回R<；ν₀+X/ν₀

可以这样编写Numpy格式副本：

def getMySamplesXpRj(rawRanCount, mu, sigma, 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，该算法的拒绝率非常低

此外，如果你愿意考虑到一些复杂的情况，也有一些关于被截断的高斯分布的Ziggurat algorithm的文献，例如尼古拉斯·肖邦在ENSAE-CREST的2012 arXiv 1201.6140 paper

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

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