python之如何以有效的方式截断 numpy/scipy 指数分布

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

我的想法是从截断为 1(下限)和 10(上限)的指数分布中抽取随机样本。我想要由此产生的有界指数分布。期望平均值为 2.5。我怎样才能有效地做到这一点？

请您参考如下方法:

有两种方法:

首先是生成一个指数分布的随机变量，然后将值限制在(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