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

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

第一种是生成指数分布的随机变量，然后将值限制为(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 ，则需要求解得到预期均值的scale参数 .

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