code: def probability_of_n_purchases_up_to_time(self, t, n):
          r
          r, alpha, a, b = self._unload_params("r", "alpha", "a", "b")
  
          first_term = (
              beta(a, b + n)
              / beta(a, b)
              * gamma(r + n)
              / gamma(r)
              / gamma(n + 1)
              * (alpha / (alpha + t)) ** r
              * (t / (alpha + t)) ** n
          )
  
          if n > 0:
              j = np.arange(0, n)
              finite_sum = (gamma(r + j) / gamma(r) / gamma(j + 1) * (t / (alpha + t)) ** j).sum()
              second_term = beta(a + 1, b + n - 1) / beta(a, b) * (1 - (alpha / (alpha + t)) ** r * finite_sum)
          else:
              second_term = 0
          return first_term + second_term
docs: r"""
          Compute the probability of n purchases.
  
           .. math::  P( N(t) = n | \text{model} )
  
          where N(t) is the number of repeat purchases a customer makes in t
          units of time.
  
          Parameters
          ----------
          t: float
              number units of time
          n: int
              number of purchases
  
          Returns
          -------
          float:
              Probability to have n purchases up to t units of time