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

1 回答

• 10

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

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

  

  回复于 2024-05-02T07:35:32+08:00