diff --git a/yli/distributions.py b/yli/distributions.py
index a8daf93..981aaa7 100644
--- a/yli/distributions.py
+++ b/yli/distributions.py
@@ -14,7 +14,6 @@
 #   You should have received a copy of the GNU Affero General Public License
 #   along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
-import mpmath
 import numpy as np
 from scipy import integrate, optimize, stats
 
@@ -50,16 +49,18 @@ class betarat_gen(stats.rv_continuous):
 	def _pdf_one(self, w, a1, b1, a2, b2):
 		"""PDF for the distribution, given by Pham-Gia"""
 		
+		from mpmath import beta, hyp2f1, power
+		
 		if w <= 0:
 			return 0
 		elif w < 1:
-			term1 = mpmath.beta(a1 + a2, b2) / (mpmath.beta(a1, b1) * mpmath.beta(a2, b2))
-			term2 = mpmath.power(w, a1 - 1)
-			term3 = mpmath.hyp2f1(a1 + a2, 1 - b1, a1 + a2 + b2, w)
+			term1 = beta(a1 + a2, b2) / (beta(a1, b1) * beta(a2, b2))
+			term2 = power(w, a1 - 1)
+			term3 = hyp2f1(a1 + a2, 1 - b1, a1 + a2 + b2, w)
 		else:
-			term1 = mpmath.beta(a1 + a2, b1) / (mpmath.beta(a1, b1) * mpmath.beta(a2, b2))
-			term2 = 1 / mpmath.power(w, a2 + 1)
-			term3 = mpmath.hyp2f1(a1 + a2, 1 - b2, a1 + a2 + b1, 1/w)
+			term1 = beta(a1 + a2, b1) / (beta(a1, b1) * beta(a2, b2))
+			term2 = 1 / power(w, a2 + 1)
+			term3 = hyp2f1(a1 + a2, 1 - b2, a1 + a2 + b1, 1/w)
 		
 		return float(term1 * term2 * term3)
 	
@@ -69,17 +70,19 @@ class betarat_gen(stats.rv_continuous):
 	def _cdf_one(self, w, a1, b1, a2, b2):
 		"""PDF for the distribution, given by Pham-Gia"""
 		
+		from mpmath import beta, hyper, power
+		
 		if w <= 0:
 			return 0
 		elif w < 1:
-			term1 = mpmath.beta(a1 + a2, b2) / (mpmath.beta(a1, b1) * mpmath.beta(a2, b2))
-			term2 = mpmath.power(w, a1) / a1
-			term3 = mpmath.hyper([a1, a1 + a2, 1 - b1], [a1 + 1, a1 + a2 + b2], w)
+			term1 = beta(a1 + a2, b2) / (beta(a1, b1) * beta(a2, b2))
+			term2 = power(w, a1) / a1
+			term3 = hyper([a1, a1 + a2, 1 - b1], [a1 + 1, a1 + a2 + b2], w)
 			return float(term1 * term2 * term3)
 		else:
-			term1 = mpmath.beta(a1 + a2, b1) / (mpmath.beta(a1, b1) * mpmath.beta(a2, b2))
-			term2 = 1 / (a2 * mpmath.power(w, a2))
-			term3 = mpmath.hyper([a2, a1 + a2, 1 - b2], [a2 + 1, a1 + a2 + b1], 1/w)
+			term1 = beta(a1 + a2, b1) / (beta(a1, b1) * beta(a2, b2))
+			term2 = 1 / (a2 * power(w, a2))
+			term3 = hyper([a2, a1 + a2, 1 - b2], [a2 + 1, a1 + a2 + b1], 1/w)
 			return 1 - float(term1 * term2 * term3)
 	
 	def _cdf(self, w, a1, b1, a2, b2):
@@ -88,8 +91,10 @@ class betarat_gen(stats.rv_continuous):
 	def _munp_one(self, k, a1, b1, a2, b2):
 		"""Moments of the distribution, given by Weekend Editor"""
 		
-		term1 = mpmath.rf(a1, k) * mpmath.rf(a2 + b2 - k, k)
-		term2 = mpmath.rf(a1 + b1, k) * mpmath.rf(a2 - k, k)
+		from mpmath import rf
+		
+		term1 = rf(a1, k) * rf(a2 + b2 - k, k)
+		term2 = rf(a1 + b1, k) * rf(a2 - k, k)
 		return float(term1 / term2)
 	
 	def _munp(self, k, a1, b1, a2, b2):
@@ -128,16 +133,18 @@ class betaoddsrat_gen(stats.rv_continuous):
 	def _pdf_one(self, w, a1, b1, a2, b2):
 		"""PDF for the distribution, given by Hora & Kelley"""
 		
+		from mpmath import beta, hyp2f1, power
+		
 		if w <= 0:
 			return 0
 		elif w < 1:
-			term1 = mpmath.beta(a1 + a2, b1 + b2) / (mpmath.beta(a1, b1) * mpmath.beta(a2, b2))
-			term2 = mpmath.power(w, a1 - 1)
-			term3 = mpmath.hyp2f1(a1 + b1, a1 + a2, a1 + a2 + b1 + b2, 1 - w)
+			term1 = beta(a1 + a2, b1 + b2) / (beta(a1, b1) * beta(a2, b2))
+			term2 = power(w, a1 - 1)
+			term3 = hyp2f1(a1 + b1, a1 + a2, a1 + a2 + b1 + b2, 1 - w)
 		else:
-			term1 = mpmath.beta(a1 + a2, b1 + b2) / (mpmath.beta(a1, b1) * mpmath.beta(a2, b2))
-			term2 = mpmath.power(w, -a2 - 1)
-			term3 = mpmath.hyp2f1(a2 + b2, a1 + a2, a1 + a2 + b1 + b2, 1 - mpmath.power(w, -1))
+			term1 = beta(a1 + a2, b1 + b2) / (beta(a1, b1) * beta(a2, b2))
+			term2 = power(w, -a2 - 1)
+			term3 = hyp2f1(a2 + b2, a1 + a2, a1 + a2 + b1 + b2, 1 - power(w, -1))
 		
 		return float(term1 * term2 * term3)
 	
@@ -147,26 +154,28 @@ class betaoddsrat_gen(stats.rv_continuous):
 	def _cdf_one_infsum(self, w, a1, b1, a2, b2):
 		"""CDF for the distribution, by truncating infinite sum given by Hora & Kelly"""
 		
+		from mpmath import beta, betainc, factorial, power, rf
+		
 		if w <= 0:
 			return 0
 		elif w < 1:
-			k = mpmath.beta(a1 + a2, b1 + b2) / (mpmath.beta(a1, b1) * mpmath.beta(a2, b2))
+			k = beta(a1 + a2, b1 + b2) / (beta(a1, b1) * beta(a2, b2))
 
 			inf_sum = 0
 			for j in range(0, BETA_ODDSRATIO_CDF_TERMS[0]):
-				kj1 = mpmath.rf(a1 + b1, j) * mpmath.rf(a1 + a2, j)
-				kj2 = mpmath.rf(a1 + a2 + b1 + b2, j) * mpmath.factorial(j)
-				inf_sum += kj1/kj2 * mpmath.betainc(a1, j + 1, 0, w)
+				kj1 = rf(a1 + b1, j) * rf(a1 + a2, j)
+				kj2 = rf(a1 + a2 + b1 + b2, j) * factorial(j)
+				inf_sum += kj1/kj2 * betainc(a1, j + 1, 0, w)
 			
 			return k * inf_sum
 		else:
-			k = mpmath.beta(a1 + a2, b1 + b2) / (mpmath.beta(a1, b1) * mpmath.beta(a2, b2))
+			k = beta(a1 + a2, b1 + b2) / (beta(a1, b1) * beta(a2, b2))
 			
 			inf_sum = 0
 			for j in range(0, BETA_ODDSRATIO_CDF_TERMS[0]):
-				kj1 = mpmath.rf(a2 + b2, j) * mpmath.rf(a1 + a2, j)
-				kj2 = mpmath.rf(a1 + a2 + b1 + b2, j) * mpmath.factorial(j)
-				inf_sum += kj1/kj2 * mpmath.betainc(a2, j + 1, 0, mpmath.power(w, -1))
+				kj1 = rf(a2 + b2, j) * rf(a1 + a2, j)
+				kj2 = rf(a1 + a2 + b1 + b2, j) * factorial(j)
+				inf_sum += kj1/kj2 * betainc(a2, j + 1, 0, power(w, -1))
 			
 			return 1 - k * inf_sum
 	
@@ -235,8 +244,10 @@ class betaoddsrat_gen(stats.rv_continuous):
 	def _munp_one(self, k, a1, b1, a2, b2):
 		"""Moments of the distribution, given by Hora & Kelley"""
 		
-		term1 = mpmath.rf(a1, k) * mpmath.rf(b2, k)
-		term2 = mpmath.rf(b1 - k, k) * mpmath.rf(a2 - k, k)
+		from mpmath import rf
+		
+		term1 = rf(a1, k) * rf(b2, k)
+		term2 = rf(b1 - k, k) * rf(a2 - k, k)
 		return float(term1 / term2)
 	
 	def _munp(self, k, a1, b1, a2, b2):