title: numpy to generate discrete probability distribution
body: <p>I'm following a code example I found at <a href="http://docs.scipy.org/doc/scipy/reference/tutorial/stats.html#subclassing-rv-discrete" rel="nofollow noreferrer">http://docs.scipy.org/doc/scipy/reference/tutorial/stats.html#subclassing-rv-discrete</a> for implementing a random number generator for discrete values of a normal distribution. The exact example (not surprisingly) works quite well, but if I modify it to allow only left or right-tailed results, the distribution around 0 should is too low (bin zero should contain more values). I must have hit a boundary condition, but am unable to work it out. Am I missing something?</p>
  
  <p>This is the result of counting the random numbers per bin:</p>
  
  <pre><code>np.bincount(rvs) [1082 2069 1833 1533 1199  837  644  376  218  111   55   20   12    7    2 2]
  </code></pre>
  
  <p>This is the histogram:</p>
  
  <p><img src="https://i.stack.imgur.com/Vl2Dw.png" alt="enter image description here"></p>
  
  <pre><code>from scipy import stats
  
  np.random.seed(42)
  
  def draw_discrete_gaussian(rng, tail='both'):
      # number of integer support points of the distribution minus 1
      npoints = rng if tail == 'both' else rng * 2
      npointsh = npoints / 2
      npointsf = float(npoints)
      # bounds for the truncated normal
      nbound = 4
      # actual bounds of truncated normal
      normbound = (1+1/npointsf) * nbound
      # integer grid
      grid = np.arange(-npointsh, npointsh+2, 1)
      # bin limits for the truncnorm
      gridlimitsnorm = (grid-0.5) / npointsh * nbound
      # used later in the analysis
      gridlimits = grid - 0.5
      grid = grid[:-1]
      probs = np.diff(stats.truncnorm.cdf(gridlimitsnorm, -normbound, normbound))
      gridint = grid
  
      normdiscrete = stats.rv_discrete(values=(gridint, np.round(probs, decimals=7)), name='normdiscrete')
      # print 'mean = %6.4f, variance = %6.4f, skew = %6.4f, kurtosis = %6.4f'% normdiscrete.stats(moments =  'mvsk')
      rnd_val = normdiscrete.rvs()
      if tail == 'both':
          return rnd_val
      if tail == 'left':
          return -abs(rnd_val)
      elif tail == 'right':
          return abs(rnd_val)
  
  
  rng = 15
  tail = 'right'
  rvs = [draw_discrete_gaussian(rng, tail=tail) for i in xrange(10000)]
  
  if tail == 'both':
      rng_min = rng / -2.0
      rng_max = rng / 2.0
  elif tail == 'left':
      rng_min = -rng
      rng_max = 0
  elif tail == 'right':
      rng_min = 0
      rng_max = rng
  
  gridlimits = np.arange(rng_min-.5, rng_max+1.5, 1)
  print gridlimits
  f, l = np.histogram(rvs, bins=gridlimits)
  
  # cheap way of creating histogram
  import matplotlib.pyplot as plt
  %matplotlib inline
  
  bins, edges = f, l
  left,right = edges[:-1],edges[1:]
  X = np.array([left, right]).T.flatten()
  Y = np.array([bins, bins]).T.flatten()
  
  # print 'rvs', rvs
  print 'np.bincount(rvs)', np.bincount(rvs)
  
  plt.plot(X,Y)
  plt.show()
  </code></pre>