numpy 生成离散概率分布答案

【问题标题】：numpy to generate discrete probability distributionnumpy 生成离散概率分布
【发布时间】：2014-03-29 16:05:20
【问题描述】：

我正在关注我在http://docs.scipy.org/doc/scipy/reference/tutorial/stats.html#subclassing-rv-discrete 找到的代码示例，用于为正态分布的离散值实现随机数生成器。确切的例子（不足为奇）效果很好，但如果我修改它以只允许左尾或右尾结果，0 附近的分布应该太低（bin 0 应该包含更多值）。我一定遇到了边界条件，但无法解决。我错过了什么吗？

这是对每个 bin 的随机数进行计数的结果：

np.bincount(rvs) [1082 2069 1833 1533 1199  837  644  376  218  111   55   20   12    7    2 2]

这是直方图：

from scipy import stats

np.random.seed(42)

def draw_discrete_gaussian(rng, tail='both'):
    # number of integer support points of the distribution minus 1
    npoints = rng if tail == 'both' else rng * 2
    npointsh = npoints / 2
    npointsf = float(npoints)
    # bounds for the truncated normal
    nbound = 4
    # actual bounds of truncated normal
    normbound = (1+1/npointsf) * nbound
    # integer grid
    grid = np.arange(-npointsh, npointsh+2, 1)
    # bin limits for the truncnorm
    gridlimitsnorm = (grid-0.5) / npointsh * nbound
    # used later in the analysis
    gridlimits = grid - 0.5
    grid = grid[:-1]
    probs = np.diff(stats.truncnorm.cdf(gridlimitsnorm, -normbound, normbound))
    gridint = grid

    normdiscrete = stats.rv_discrete(values=(gridint, np.round(probs, decimals=7)), name='normdiscrete')
    # print 'mean = %6.4f, variance = %6.4f, skew = %6.4f, kurtosis = %6.4f'% normdiscrete.stats(moments =  'mvsk')
    rnd_val = normdiscrete.rvs()
    if tail == 'both':
        return rnd_val
    if tail == 'left':
        return -abs(rnd_val)
    elif tail == 'right':
        return abs(rnd_val)


rng = 15
tail = 'right'
rvs = [draw_discrete_gaussian(rng, tail=tail) for i in xrange(10000)]

if tail == 'both':
    rng_min = rng / -2.0
    rng_max = rng / 2.0
elif tail == 'left':
    rng_min = -rng
    rng_max = 0
elif tail == 'right':
    rng_min = 0
    rng_max = rng

gridlimits = np.arange(rng_min-.5, rng_max+1.5, 1)
print gridlimits
f, l = np.histogram(rvs, bins=gridlimits)

# cheap way of creating histogram
import matplotlib.pyplot as plt
%matplotlib inline

bins, edges = f, l
left,right = edges[:-1],edges[1:]
X = np.array([left, right]).T.flatten()
Y = np.array([bins, bins]).T.flatten()

# print 'rvs', rvs
print 'np.bincount(rvs)', np.bincount(rvs)

plt.plot(X,Y)
plt.show()

【问题讨论】：

查看图表，在我看来 bin 0 包含从 -0.5 到 0.5 的所有内容。如果是这样，它大约是下一个垃圾箱的一半也就不足为奇了。您没有从该 bin 的左半部分生成结果。
@user2357112：我可能是错的，但我认为这只是由于可视化（它以 bin 编号为中心，而实际上 bin 以 +0.5 为界）。如果我做gridlimits = np.arange(rng_min, rng_max+2, 1)，这是同一张图。
我也认为@user235711 是对的。当您采用 abs 时，您正在组合 probs 的负箱和正箱。检查从零开始的 bin 的长度是否与其他 bin 的组合长度相同。我只会为右或左取正确的截断法线，即从零开始或结束。
你所说的两个听起来都令人信服（尽管我认为我通过乘以支持点的数量npoints = rng if tail == 'both' else rng * 2 绕过了这个半拆分箱 0）...
感谢您的两位 cmets，我可能已经找到了解决方案 - 感谢您的帮助 (+1x2)。

标签： python numpy scipy

【解决方案1】：

我尝试根据@user333700 和@user235711 的cmets 回答我自己的问题：

我在normdiscrete = ...之前插入方法

if tail == 'right':
    gridint = gridint[npointsh:]
    probs = probs[npointsh:]
    s = probs.sum()
    probs = probs / s
elif tail == 'left':
    gridint = gridint[0: npointsh]
    probs = probs[0: npointsh]
    s = probs.sum()
    probs = probs / s

生成的直方图和看起来更好：

【讨论】：