Noticed an error in line 59:
mom.append(np.sum(xpf)) # compute moment: sum(x^p * f(x))
moms is defined before as moms = [], therefore "mom" in line 59 should be changed to "moms".
I'm also having trouble understanding the use for randvar, if you could kindly explain. Thank you.