def sample(self, n):
lo, hi = self.bounds<EOL>C = (<NUM_LIT:1> + self.alpha)/(hi**(<NUM_LIT:1> + self.alpha) - lo**(<NUM_LIT:1> + self.alpha))<EOL>u = np.random.random(n)<EOL>a = self.alpha<EOL>return ((a+<NUM_LIT:1>) * (u/C + (lo**(a+<NUM_LIT:1>) / (a+<NUM_LIT:1>))))**(<NUM_LIT:1>/(a+<NUM_LIT:1>))<EOL>
cdf = C * int(x**a)|x=lo..x

    = C * [x**(a+1) / (a+1)] | x=lo..x
    = C * ([x**(a+1) / (a+1)] - [lo**(a+1) / (a+1)])
 u  =

 u/C + [lo**(a+1) / (a+1)] = x**(a+1) / (a+1)
 (a+1) * (u/C + [lo**(a+1) / (a+1)]) = x**(a+1)
 [(a+1) * (u/C + [lo**(a+1) / (a+1)])] ** (1/(a+1)) = x