使用 numpy.random.normal 时如何指定上限和下限

Posted 2023-02-23

【中文标题】使用 numpy.random.normal 时如何指定上限和下限

【英文标题】：How to specify upper and lower limits when using numpy.random.normal

【发布时间】：2013-08-28 19:18:01

【问题描述】：

我希望能够从仅介于 0 和 1 之间的正态分布中选择值。在某些情况下，我希望能够基本上只返回一个完全随机的分布，而在其他情况下，我希望返回值呈高斯形状。

目前我正在使用以下功能：

def blockedgauss(mu,sigma):
    while True:
        numb = random.gauss(mu,sigma)
        if (numb > 0 and numb < 1):
            break
    return numb

它从正态分布中选择一个值，然后在它超出 0 到 1 的范围时将其丢弃，但我觉得必须有更好的方法来做到这一点。

【问题讨论】：

如果你“阻止”值 1，它仍然是高斯分布吗？

它不会是高斯分布，但在某些情况下我不想要高斯分布。我想返回一个可在随机分布（从非常宽的高斯中挑选）到非常接近 delta 函数（高斯变得非常窄）之间可调的分布

【参考方案1】：

听起来您想要truncated normal distribution。 使用 scipy，您可以使用 scipy.stats.truncnorm 从这样的分布中生成随机变量：

import matplotlib.pyplot as plt
import scipy.stats as stats

lower, upper = 3.5, 6
mu, sigma = 5, 0.7
X = stats.truncnorm(
    (lower - mu) / sigma, (upper - mu) / sigma, loc=mu, scale=sigma)
N = stats.norm(loc=mu, scale=sigma)

fig, ax = plt.subplots(2, sharex=True)
ax[0].hist(X.rvs(10000), normed=True)
ax[1].hist(N.rvs(10000), normed=True)
plt.show()

上图为截断正态分布，下图为均值相同mu、标准差相同sigma的正态分布。

【讨论】：

为了方便，我使用 sigma = (upper-lower)/sigma_span。对于 sigma_span = 4（+2s 和 -2s），正态分布的 95% 的值将落入界限。 sigma_span=1 产生非常广泛的分布，较高的值 (>6) 与正态分布几乎无法区分（大多数值不会达到边界）。

【参考方案2】：

我在寻找一种方法来返回从 0 到 1 之间截断的正态分布中采样的一系列值（即概率）时遇到了这篇文章。为了帮助遇到同样问题的其他人，我只想指出 scipy.stats.truncnorm 具有内置功能“.rvs”。

所以，如果您想要 100,000 个样本，平均值为 0.5，标准差为 0.1：

import scipy.stats
lower = 0
upper = 1
mu = 0.5
sigma = 0.1
N = 100000

samples = scipy.stats.truncnorm.rvs(
          (lower-mu)/sigma,(upper-mu)/sigma,loc=mu,scale=sigma,size=N)

这给出了与 numpy.random.normal 非常相似的行为，但在所需的范围内。使用内置函数将比循环收集样本快得多，尤其是对于较大的 N 值。

【讨论】：

当下限为负时，这似乎永远不会返回负数。

scipy 是否支持其他分发类型？

@tagoma 如果你提高 sigma，你应该会看到负数。

【参考方案3】：

如果有人想要仅使用 numpy 的解决方案，这里有一个使用 normal 函数和 clip 的简单实现（MacGyver 的方法）：

    import numpy as np
    def truncated_normal(mean, stddev, minval, maxval):
        return np.clip(np.random.normal(mean, stddev), minval, maxval)

编辑：不要使用这个！！这是你不应该这样做的方式！！例如，a = truncated_normal(np.zeros(10000), 1, -10, 10) 可能看起来可行，但是b = truncated_normal(np.zeros(10000), 100, -1, 1) 将绝对不会绘制截断的法线，您可以在下面的直方图中看到：

对不起，希望没有人受伤！我想教训是，不要试图在编码方面模仿 MacGyver ...... 干杯， 安德烈斯

【讨论】：

麦盖弗是谁？

这是一个大胆的角色，可以用普通的常规项目构建复杂的系统来解决困难的情况。这个答案有点模因，请尽情参考 MacGyver！

【参考方案4】：

我通过以下方式制作了一个示例脚本。它展示了如何使用 API 来实现我们想要的功能，例如生成具有已知参数的样本，如何计算 CDF、PDF 等。我还附上一张图片来展示这一点。

#load libraries   
import scipy.stats as stats

#lower, upper, mu, and sigma are four parameters
lower, upper = 0.5, 1
mu, sigma = 0.6, 0.1

#instantiate an object X using the above four parameters,
X = stats.truncnorm((lower - mu) / sigma, (upper - mu) / sigma, loc=mu, scale=sigma)

#generate 1000 sample data
samples = X.rvs(1000)

#compute the PDF of the sample data
pdf_probs = stats.truncnorm.pdf(samples, (lower-mu)/sigma, (upper-mu)/sigma, mu, sigma)

#compute the CDF of the sample data
cdf_probs = stas.truncnorm.cdf(samples, (lower-mu)/sigma, (upper-mu)/sigma, mu, sigma)

#make a histogram for the samples
plt.hist(samples, bins= 50,normed=True,alpha=0.3,label='histogram');

#plot the PDF curves 
plt.plot(samples[samples.argsort()],pdf_probs[samples.argsort()],linewidth=2.3,label='PDF curve')

#plot CDF curve        
plt.plot(samples[samples.argsort()],cdf_probs[samples.argsort()],linewidth=2.3,label='CDF curve')


#legend
plt.legend(loc='best')

【参考方案5】：

实际上，您可以对数据进行规范化，然后将其传输到您需要的范围。不好意思第一次用，不知道怎么直接显示图片 the function is shown

【参考方案6】：

我已经使用 numpy 测试了一些解决方案。通过反复试验，我发现± variation 除以 3 是标准差的一个很好的猜测。

以下是一些示例：

基础知识

import numpy as np
import matplotlib.pyplot as plt

val_min = 1000
val_max = 2000
variation = (val_max - val_min)/2
std_dev = variation/3
mean = (val_max + val_min)/2
dist_normal = np.random.normal(mean, std_dev,  1000)
print('Normal distribution\n\tMin: 0:.2f, Max: 1:.2f'
      .format(dist_normal.min(), dist_normal.max()))
plt.hist(dist_normal, bins=30)
plt.show()

比较案例

import numpy as np
import matplotlib.pyplot as plt

val_min = 1400
val_max = 2800
variation = (val_max - val_min)/2
std_dev = variation/3
mean = (val_max + val_min)/2
fig, ax = plt.subplots(3, 3)
plt.suptitle("Histogram examples by Davidson Lima (github.com/davidsonlima)", 
             fontweight='bold')
i = 0
j = 0
pos = 1
while (i < 3):
    while (j < 3):
        dist_normal = np.random.normal(mean, std_dev,  1000)
        max_min = 'Min: 0:.2f, Max: 1:.2f'.format(dist_normal.min(), dist_normal.max())
        ax[i, j].hist(dist_normal, bins=30, label='Dist' + str(pos))
        ax[i, j].set_title('Normal distribution ' + str(pos))
        ax[i, j].legend()
        ax[i, j].text(mean, 0, max_min, horizontalalignment='center', color='white',
                      bbox='facecolor': 'red', 'alpha': 0.5)
        print('Normal distribution 0\n\tMin: 1:.2f, Max: 2:.2f'
              .format(pos, dist_normal.min(), dist_normal.max()))
        j += 1
        pos += 1
    j = 0
    i += 1
plt.show()

如果有人对 numpy 有更好的方法，请在下面发表评论。