<p>最简单的方法可能是利用Scipy提供的截断正态分布</p>
<p>这给出了以下代码,其中ν(nu)作为标准高斯分布的变量,τ(tau)映射到该分布上的ν<sub>0</sub>。此函数返回包含ranCount lognormal变量的Numpy数组:</p>
<pre><code>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
</code></pre>
<p>如果出于某种原因,这是不合适的,那么一些常用于高斯变量的技巧,例如<a href="https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform" rel="nofollow noreferrer">Box-Muller</a>对截断分布不起作用,但是我们可以始终求助于一个一般原则:<a href="https://en.wikipedia.org/wiki/Inverse_transform_sampling" rel="nofollow noreferrer">Inverse Transform Sampling</a>定理</p>
<p>所以我们通过变换均匀变量来生成变量的累积概率。我们信任Scipy,使用它的<em>erf</em>误差函数的逆函数,从概率返回到x空间值</p>
<p>这提供了类似于以下Python代码的内容(无需任何优化尝试):</p>
<pre><code>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
</code></pre>
<h2 id="alternatives">备选方案:</h2>
<p>我们可以利用大量与<a href="https://en.wikipedia.org/wiki/Truncated_normal_distribution" rel="nofollow noreferrer">Truncated Gaussian distribution</a>相关的材料</p>
<p>Zdravko Botev和Pierre L'Ecuyer最近(2016年)就这一主题发表了一篇文章。本文提供了一个指向公共可用的<a href="https://cran.r-project.org/web/packages/TruncatedNormal" rel="nofollow noreferrer">R source code</a>的指针。有些材料非常陈旧,例如1986年吕克·德夫罗耶的书:<a href="http://www.eirene.de/Devroye.pdf" rel="nofollow noreferrer">Non-Uniform Random Variate Generation</a></p>
<p>例如,一种可能的基于拒绝的方法:如果τ(tau)映射到标准高斯曲线上的ν<sub>0</sub>,则单位高斯分布类似于exp(-ν<sup>2</sup>/2)。如果我们写ν=ν<sub>0</sub>+δ,这与exp(-δ<sup>2</sup>/2)*exp(-ν<sub>0</sub>*δ)成正比</p>
<p>其思想是通过参数ν<sub>0</sub>的<strong>指数分布</strong>来近似超过ν<sub>0</sub>的精确分布。请注意,精确分布始终低于近似分布。然后我们可以随机接受相对便宜的指数变量,概率为exp(-δ<sup>2</sup>/2)</p>
<p>我们可以在文献中选择一个等价的算法。在Devroye书第九章第382页中,有一些伪代码:</p>
<p>重复
生成独立的指数随机变量X和Y
直到X<sup>2</sup><;=2*ν<sub>0</sub><sup>2</sup>*Y</p>
<p>返回R<;ν<sub>0</sub>+X/ν<sub>0</sub></p>
<p>可以这样编写Numpy格式副本:</p>
<pre><code>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
</code></pre>
<p>根据Botev-L'Ecuyer论文中的表3,该算法的拒绝率非常低</p>
<p>此外,如果你愿意考虑到一些复杂的情况,也有一些关于被截断的高斯分布的<a href="https://en.wikipedia.org/wiki/Ziggurat_algorithm" rel="nofollow noreferrer"><em>Ziggurat</em> algorithm</a>的文献,例如尼古拉斯·肖邦在ENSAE-CREST的2012 <a href="https://arxiv.org/pdf/1201.6140" rel="nofollow noreferrer">arXiv 1201.6140 paper</a></p>
<p><strong>旁注:</strong>对于Python的最新版本,似乎可以直接使用希腊字母作为变量名,σ代替σ,τ代替τ,就像在统计书籍中一样:</p>
<pre><code>$ python3
Python 3.9.6 (default, Jun 29 2021, 00:00:00)
>>>
>>> σ = 2
>>> τ = 7
>>>
>>> στ = σ * τ
>>>
>>> στ + 1
15
>>>
</code></pre>