From e1fd73e095eac67b7806aab213f72127f94c3447 Mon Sep 17 00:00:00 2001
From: RunasSudo <runassudo@yingtongli.me>
Date: Thu, 6 Oct 2022 18:30:25 +1100
Subject: [PATCH] Add test case for beta_ratio.pdf/cdf/mean

---
 .gitignore                      |   2 +
 tests/__init__.py               |   0
 tests/beta_ratio_vs_jsaffer.npy | Bin 0 -> 1728 bytes
 tests/test_beta_ratio.py        |  80 ++++++++++++++++++++++++++++++++
 4 files changed, 82 insertions(+)
 create mode 100644 .gitignore
 create mode 100644 tests/__init__.py
 create mode 100644 tests/beta_ratio_vs_jsaffer.npy
 create mode 100644 tests/test_beta_ratio.py

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 0000000000000000000000000000000000000000..71b4851e89bb8384ef83705a20aa83f35c57b1eb
GIT binary patch
literal 1728
zcmbW0`#Ti)9>-h7E=M*}?y|I5vO>xh&fB3la>?VMs1%u^tYOBI+pgSFO(VCIETtHG
za*Y}pZ8@Y^L>e;Laang}zA-L?Y-i5%obwNy_b;#K{r=_ie4h97&UE#3ar6FcUBtQ=
z{QwI6G);fMj=t5I@AQpy^v_Uf;WWRHlT=y&<v+Q-Uof2_$?1W9p%jV##>_~^)WpQV
zNGDF`|CHRi|0>gXkWnv4M=+bl&U##l@?(eM|4e9rfg%6wi={RQ1U#u4SuQAF*xCi?
z_M=JB){5CU1a6BUZ8?4vadsvAY_my}ZZO+w^I`_ZJSl#j@&|BtP=1;&nZ+KZMAl7p
z5qQZ|6GE9s?vi)CXX*lYD#kS9_(jxR4sNJ*UP4T}uJ(dp3E{~TMOOmEXtsWpUtT4~
zUk8?-Z|oK0#)5|xyGM+W*;D6atHp3+PHNo=5#vqLTKoZdF`idh)A-CK6!iJKpIu%=
zi|L>l-*piKJ*i&S84Hm1)?Ns0o`-!(x!GZ%2#&XlEz{*i;7BDzwCl}*{o4wMws#ia
zf4ZB}Wd8xZAFpnyEEB?4L;djft{E&udK_?Gox=9CufH?gHHk!T^K-|k;~12CoTdHi
zD2$tmlAX0ia4*&P5rsL7>E_LxbfqB#U2Bnb@*KeX8_W~#<^6Dpa%!kI;UUdYU60w{
zi@jwEJG`s7Xi<(~xRiIpB<182`!`)U&d-(SX>=kjdOfE_sU1$&-&gE5YX$E_%X*EB
zmq_4uoJ%x*fhUEDFHf#-M0lUO+<o~5SO|~YE@Au%fvDiLt~Uq8ycI*c$Xc*Hw{S!Q
zkD&XnWR(O)HAatB^H>p8h_XqmP&chWQ6%%5tiA`ZF#r3(=Kd1UYDcTnOp75ZjH<|e
zeGj4}eRubYyU6(HS6G%=fVe%=w!@}*XeSF@^2WCiQeL`M<;NQcoNwGZr<w)&)%JHk
zE~X>JNaN&%FRvputo6^QbqpM!ulH=NxeO1tcOD$6OSluRR_t{<2?MIqdjld8pm51s
zEt3(C%uOM)4drpT5$jU-d0{lv9K8HWCnE8W{Z)+1Q9oegxNmHnGab5lv*xZoVc^G;
zsBQbGNFLZ_r>+qUvz?~mvu;7)b@YU6=sW|u|B0IE^Z=Y;=UA^X{ZXJ-<~IEN6m0yJ
z@`ZN3*cl|$qWT=ePuwq;z3PwRfBU%mes|ql`ZjRqj?%C4T%fXhMk{6PECwSR!W{Y0
z(7ot(Xva`0uJ!IqFg%ro?t{N1beH5qepUFnYFQDya%^^t?JI+WSkZHG`5}%**XC=D
zvl0HlM0)<%Q<SP{AN^F^0FwH__R+g0)Tic~R$9G6hoJA(w{5?n>AbI^W&CT5s0$XR
zO1hx*PsT>+pSv+Q|CO5bEiTT7ICvM{>cy4j;?BC;eQ2}_3-QT&4_UEG-SeBjLyp1F
zr7#CzTdOJ)_jC{$qj`Q7<wKBZ(lWnrb6C<3Yc)LP4?G*c&~VXf1R5rLeK#45g5@g;
zt(Y4{j6F+Pwt5Wa*C&04g2o}W$LXQ6$^=SfRs`8~6Zrj8V`8w!B>G3^Vjqr7VmhP$
zhm7-6uvZAoIXF3ms(_u@ioVl`mNI;7-#U#-tpb6y@eC}l-TUIml^LA37@l9`&meWv
zoL-8N0352D0+lWR+tuftW{m)@Ju(_ALIE@bFLzP53*oIt2_EqfVzzQXmJ=_8P`tg{
zy+nwiimq6hHX-I;e>S;uMu-|l=)`7O0#+zTMNW_4O-C%(V?V*r@DrCkhY4~N6m!P?
z2zZ?a{6ZRmhpcAat&0TNjU8{7nFJlrSU!SW0t?oGExAPmZ<UTX>Q)eBWUCI|XA}6d
zbjqxs5p1Zhs*-Ld82I|_@y%@n@|kTfqq_(?!+sIE^bp)PcsY{CBe=tLkdELJe5;wh
zZET2Ob$YHqJVMazR8+$rCz2hxQ;ntwK31Rm!*@n<_p9)~c?k*nUG5Rp4+Jq*;{~O&
z1VI%KPHz<vtfePC^_eHo=Pp#nEfADyUFu9&BnbS=;X~gq5!ktEiy)T#X-jomldw3-
z@-Azcz%Q6PYO_M{^ZlW*RtY-<mn<Dt37ENau51Z&3U*O6KN5r+B^N|W=$GWUu2sS}
Ry3QfWYXn-Q{p=$W{u}kDQ-c5i

literal 0
HcmV?d00001

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