diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..8d35cb3
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,2 @@
+__pycache__
+*.pyc
diff --git a/tests/__init__.py b/tests/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/tests/beta_ratio_vs_jsaffer.npy b/tests/beta_ratio_vs_jsaffer.npy
new file mode 100644
index 0000000..71b4851
Binary files /dev/null and b/tests/beta_ratio_vs_jsaffer.npy differ
diff --git a/tests/test_beta_ratio.py b/tests/test_beta_ratio.py
new file mode 100644
index 0000000..605bfb4
--- /dev/null
+++ b/tests/test_beta_ratio.py
@@ -0,0 +1,80 @@
+#   scipy-yli: Helpful SciPy utilities and recipes
+#   Copyright © 2022  Lee Yingtong Li (RunasSudo)
+#
+#   This program is free software: you can redistribute it and/or modify
+#   it under the terms of the GNU Affero General Public License as published by
+#   the Free Software Foundation, either version 3 of the License, or
+#   (at your option) any later version.
+#
+#   This program is distributed in the hope that it will be useful,
+#   but WITHOUT ANY WARRANTY; without even the implied warranty of
+#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#   GNU Affero General Public License for more details.
+#
+#   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 numpy as np
+
+import yli
+
+def test_beta_ratio_vs_jsaffer_pdf():
+	"""Compare beta_ratio.pdf with result from https://github.com/jsaffer/beta_quotient_distribution"""
+	
+	# Define the example beta distribution
+	a1, b1, a2, b2 = 3, 6, 12, 7
+	
+	# Compute PDF
+	x = np.linspace(0, 2, 100)
+	dist = yli.beta_ratio(a1, b1, a2, b2)
+	y = dist.pdf(x)
+	
+	# Compare with expected values from jsaffer implementation
+	expected = np.load('beta_ratio_vs_jsaffer.npy', allow_pickle=False)[0]
+	
+	# Allow 1e-10 tolerance
+	diff = np.abs(y - expected)
+	assert (diff < 1e-10).all()
+
+def test_beta_ratio_vs_jsaffer_cdf():
+	"""Compare beta_ratio.cdf with result from https://github.com/jsaffer/beta_quotient_distribution"""
+	
+	# Define the example beta distribution
+	a1, b1, a2, b2 = 3, 6, 12, 7
+	
+	# Compute PDF
+	x = np.linspace(0, 2, 100)
+	dist = yli.beta_ratio(a1, b1, a2, b2)
+	y = dist.cdf(x)
+	
+	# Compare with expected values from jsaffer implementation
+	expected = np.load('beta_ratio_vs_jsaffer.npy', allow_pickle=False)[1]
+	
+	# Allow 1e-10 tolerance
+	diff = np.abs(y - expected)
+	assert (diff < 1e-10).all()
+
+def _gen_beta_ratio_vs_jsaffer():
+	"""Generate beta_ratio_vs_jsaffer.npy for test_beta_ratio_vs_jsaffer_pdf/cdf"""
+	
+	import beta_quotient_distribution
+	
+	a1, b1, a2, b2 = 3, 6, 12, 7
+	
+	x = np.linspace(0, 2, 100)
+	y1 = np.vectorize(lambda w: float(beta_quotient_distribution.pdf_bb_ratio(a1, a2, b1, b2, w)))(x)
+	y2 = np.vectorize(lambda w: float(beta_quotient_distribution.cdf_bb_ratio(a1, a2, b1, b2, w)))(x)
+	
+	np.save('beta_ratio_vs_jsaffer.npy', np.array([y1, y2]), allow_pickle=False)
+
+def test_beta_ratio_mean_vs_empirical():
+	"""Compare beta_ratio.mean (via beta_ratio._munp) with empirical mean"""
+	
+	# Define the example beta distribution
+	dist = yli.beta_ratio(3, 6, 12, 7)
+	
+	# Compute empirical mean
+	sample = dist.rvs(1000, random_state=31415)
+	
+	# Allow 0.01 tolerance
+	assert np.abs(dist.mean() - sample.mean()) < 0.01