如何最好地从截断正态分布中获取样本？

Question

我已经做了一些搜索，但我似乎无法找到从截断正态分布中采样的合理方法。

没有截断我正在做：

samples = [np.random.normal(loc=x,scale=d) for (x,d) in zip(X,D)]

X并D成为浮动列表。

目前我正在实施截断：

def truncnorm(loc,scale,bounds):
  s = np.random.normal(loc,scale)
  if s > bounds[1]:
    return bounds[1]
  elif s < bounds[0]:
    return bounds[0]
  return s

samples = [truncnorm(loc=x,scale=d,bounds=b) for (x,d,b) in zip(X,D,bounds)]

bounds是一个元组列表(min,max)

这种方法感觉有点尴尬，所以想知道是否有更好的方法？

Answer 1

返回边界之外的样本的值将导致太多样本落在边界上。这并不代表实际分布。边界上的值需要被拒绝并由新样本替换。这样的代码可以是：

def test_truncnorm(loc, scale, bounds):
    while True:
        s = np.random.normal(loc, scale)
        if bounds[0] <= s <= bounds[1]:
            break
    return s

由于边界很窄，这可能会非常慢。Scipy 的truncnorm可以更有效地处理此类情况。有点令人惊讶的是，边界以标准法线的函数表示，因此您的调用将是：

s = scipy.stats.truncnorm.rvs((bounds[0]-loc)/scale, (bounds[1]-loc)/scale, loc=loc, scale=scale)

请注意，当使用 numpy 的矢量化和广播时，scipy 的工作速度要快得多。一旦您习惯了这种符号，它的书写和阅读看起来也会更简单。所有样本均可一次性计算如下：

X = np.array(X)
D = np.array(D)
bounds = np.array(bounds)
samples = scipy.stats.truncnorm.rvs((bounds[:, 0] - X) / D, (bounds[:, 1] - X) / D, loc=X, scale=D)