code: def rtruncated_poisson(mu, k, size=None):
      
  
      # Calculate m
      try:
          m = max(0, np.floor(k - mu))
      except (TypeError, ValueError):
          # More than one mu
          return np.array([rtruncated_poisson(x, i)
                           for x, i in zip(mu, np.resize(k, np.size(mu)))]).squeeze()
      k -= 1
  
      # Calculate constant for acceptance probability
      C = np.exp(flib.factln(k + 1) - flib.factln(k + 1 - m))
  
      # Empty array to hold random variates
      rvs = np.empty(0, int)
  
      total_size = np.prod(size or 1)
  
      while(len(rvs) < total_size):
  
          # Propose values by sampling from untruncated Poisson with mean mu + m
          proposals = np.random.poisson(
              mu + m, (total_size * 4, np.size(m))).squeeze()
  
          # Acceptance probability
          a = C * np.array([np.exp(flib.factln(y - m) - flib.factln(y))
                            for y in proposals])
          a *= proposals > k
  
          # Uniform random variates
          u = np.random.random(total_size * 4)
  
          rvs = np.append(rvs, proposals[u < a])
  
      return np.reshape(rvs[:total_size], size)
docs: Random truncated Poisson variates with minimum value k, generated
      using rejection sampling.