在 Python 中指定对数正态分布的范围

Question

我正在使用内置函数生成随机对数正态分布。是否可以为给定的 mean 和 sigma 指定一个范围？例如，我想要 0 和 0.5 之间的随机分布 mean=0.2, sigma=0.5?

import numpy as np 

r=np.random.lognormal(mean=0.2, sigma=0.5, size=(3,3))
print([r])

Answer 1

不是直接的，但您可以生成更多点并过滤那些超出参数的点。但是，使用您提供的值，这样的点将很少：

np.random.seed(0)
z = np.random.lognormal(mean=0.2, sigma=0.5, size=10_000)

>>> np.sum(z < 0.5) / len(z)
0.0346

作为旁注，请注意 np.random.lognormal() 的参数 mean 和 sigma 指的是基础正态分布，而不是 mean 和 std 的 log-normal 点：

np.random.seed(0)
z = np.random.lognormal(mean=0.2, sigma=0.5, size=1000)
y = np.log(y)

>>> np.mean(z), np.std(z)
(1.3886515119063163, 0.7414709414626542)

>>> np.mean(y), np.std(y)
(0.2018986489272414, 0.5034218384643446)

综上所述，如果您真的想要您所要求的，您可以这样做：

shape = 3, 3
zmin, zmax = 0, 0.5

n = np.prod(shape)
zc = np.array([])
while True:
    z = np.random.lognormal(mean=0.2, sigma=0.5, size=n * 100)
    z = z[(zmin <= z) & (z < zmax)]
    zc = np.r_[zc, z]
    if len(zc) >= n:
        break
r = zc[:n].reshape(shape)

>>> r
array([[0.34078793, 0.45366409, 0.40183988],
       [0.43387773, 0.46387512, 0.30535007],
       [0.44248787, 0.32316698, 0.48600577]])

但是请注意，如果您的边界指定一个空 space.

，上面的循环（在本例中平均只完成一次）可能运行永远

编辑

为了效率，我们可以改用scipy.stats.truncnorm。不过loc, scale != 0, 1的时候好像有问题，所以这里推出自己的调整：

# convert X ~ N(0, 1) to Y ~ N(u, s):
# Y = X * s + u
# want Y in [a, b]:
# so X in [(a - u)/s, (b - u)/s]
#
# We want Z ~ logN(u, s) in [zmin, zmax]
# Z = exp(Y) and
# a, b = np.log((zmin, zmax))

def lognormal(a, b, loc=1, scale=1, size=1):
    a, b = np.log((a, b))
    ca, cb = (np.array((a, b)) - loc) / scale
    rv = truncnorm(a=ca, b=cb)
    z = np.exp(rv.rvs(size=size) * scale + loc)
    return z

示例：

zmin, zmax = 0, 0.5
z = lognormal(zmin, zmax, 0.2, 0.5, 10_000)

或者，要捕获来自 np.log(0) 的相当嘈杂的警告：

import warnings

with warnings.catch_warnings():
    warnings.filterwarnings('ignore', category=RuntimeWarning, message='divide by zero encountered in log') 
    z = lognormal(zmin, zmax, 0.2, 0.5, 10_000)

无论如何，我们现在有：

np.min(z), np.max(z)
(0.10787873738453256, 0.4999993445360746)

并且与上面的 (r) 和 shape = 10_000:

得到的分布相同

def lognormal_naive(a, b, loc=1, scale=1, size=1):
    zc = np.array([])
    while True:
        z = np.random.lognormal(mean=loc, sigma=scale, size=n * 100)
        z = z[(a <= z) & (z < b)]
        zc = np.r_[zc, z]
        if len(z) >= size:
            break
    return zc[:size]

r = lognormal_naive(zmin, zmax, 0.2, 0.5, 10_000)

fig, ax = plt.subplots(ncols=2, figsize=(7, 2))
ax[0].hist(z, bins=100)
ax[1].hist(r, bins=100);

然而，时机要好得多：

%timeit lognormal_naive(zmin, zmax, 0.2, 0.5, 10_000)
64.2 ms ± 60.3 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

%timeit lognormal(zmin, zmax, 0.2, 0.5, 10_000)
2.06 ms ± 3.61 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

在 Python 中指定对数正态分布的范围

Specifiying range for log-normal distribution in Python

python

numpy

distribution