如何以有效的方式截断 numpy/scipy 指数分布？答案

【问题标题】：How to truncate a numpy/scipy exponential distribution in an efficient way?如何以有效的方式截断 numpy/scipy 指数分布？
【发布时间】：2014-09-28 06:13:03
【问题描述】：

我目前正在构建一个神经科学实验。基本上，每 x 秒（x = inter-trial interval）会呈现 3 秒的刺激。我希望 x 相当短 (mean = 2.5) 且不可预测。

我的想法是从截断为 1（下限）和 10（上限）的指数分布中抽取随机样本。我想要得到的有界指数分布。期望均值为 2.5。我怎样才能以有效的方式做到这一点？

【问题讨论】：

如果您希望样本有界，为什么要选择指数分布？
请注意，如果从均值为 2.5 的指数分布开始，然后截断到区间 [1, 10]，则截断分布的均值不是 2.5。其实是3.25左右。

标签： python statistics scipy distribution

【解决方案1】：

有两种方法可以做到这一点：

首先是生成一个指数分布的随机变量，然后将值限制为(1,10)。

In [14]:

import matplotlib.pyplot as plt
import scipy.stats as ss
Lambda = 2.5 #expected mean of exponential distribution is lambda in Scipy's parameterization
Size = 1000
trc_ex_rv = ss.expon.rvs(scale=Lambda, size=Size)
trc_ex_rv = trc_ex_rv[(trc_ex_rv>1)&(trc_ex_rv<10)]
In [15]:

plt.hist(trc_ex_rv)
plt.xlim(0, 12)
Out[15]:
(0, 12)

In [16]:

trc_ex_rv
Out[16]:
array([...]) #a lot of numbers

当然，问题是你不会得到随机数的确切数量（这里由Size 定义）。

另一种方法是使用Inverse transform sampling，您将获得指定的准确重复次数：

In [17]:
import numpy as np
def trunc_exp_rv(low, high, scale, size):
    rnd_cdf = np.random.uniform(ss.expon.cdf(x=low, scale=scale),
                                ss.expon.cdf(x=high, scale=scale),
                                size=size)
    return ss.expon.ppf(q=rnd_cdf, scale=scale)
In [18]:

plt.hist(trunc_exp_rv(1, 10, Lambda, Size))
plt.xlim(0, 12)
Out[18]:
(0, 12)

如果您希望得到的有界分布具有给定值的预期均值，例如2.5，则需要求解产生预期均值的尺度参数。

import scipy.optimize as so
def solve_for_l(low, high, ept_mean):
    A = np.array([low, high])
    return 1/so.fmin(lambda L: ((np.diff(np.exp(-A*L)*(A*L+1)/L)/np.diff(np.exp(-A*L)))-ept_mean)**2,
                     x0=0.5,
                     full_output=False, disp=False)
def F(low, high, ept_mean, size):
    return trunc_exp_rv(low, high,
                        solve_for_l(low, high, ept_mean),
                        size)
rv_data = F(1, 10, 2.5, 1e5)
plt.hist(rv_data, bins=50)
plt.xlim(0, 12)
print rv_data.mean()

结果：

2.50386617882

【讨论】：

请注意，此分布的平均值不是 2.5，这是@user1363251 在问题第一段中要求的。
是的，我注意到了。不完全确定 OP 是否想要截断预期均值为 2.5 的指数分布，或者想要得到的有界指数分布。期望均值为 2.5。如果他能澄清一下，那是很容易做到的。期望均值可能很容易用封闭式表达。
如果您使用 hist （将输入“分组”）然后在绘图中使用 xlim，我不明白您的意思是如何截断。问题是给定一个大样本，例如1000，房车。如何过滤它们以仅获取适合给定范围的那些。房车是无序的。
>>> a = expon.rvs(size=10, scale=4)array([6.05231017, 5.71233108, 1.74126937, 0.57614087, 7.5496941 , 2.84520597, 3.06463358, 2.84783107, 1.84672096, 4.45238162])>>> np.mean(a)3.6688518790529456>>> np.sort(a)[0:5]array([0.57614087, 1.74126937, 1.84672096, 2.84520597, 2.84783107])>>> np.mean(np.sort(a)[0:5])1.9714336475094363 # different mean!

【解决方案2】：

除了@CT Zhu 的出色回答外，scipy 现在似乎内置了truncated exponential distribution。

from scipy.stats import truncexpon
r = truncexpon.rvs(b, size=1000)

【讨论】：