如何高效地截断numpy/scipy的指数分布？

Question

我现在正在做一个神经科学实验。简单来说，就是每隔一段时间（x = 试验间隔），会展示一个刺激，持续3秒。我希望这个间隔x比较短（平均值 = 2.5），而且是不可预测的。

我的想法是从一个被限制在1到10之间的指数分布中随机抽样。我希望这个限制后的指数分布的期望值能达到2.5。我该怎么高效地实现这个呢？

Answer 1

除了@CT Zhu给出的很棒的回答，看来现在scipy里已经内置了一个叫做截断指数分布的功能。

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

Answer 2

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

第一种方法是生成一个指数分布的随机数，然后把这些数限制在(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来定义）。

另一种方法是使用逆变换抽样，这样你就能得到指定数量的重复值：

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