RunasSudo/scipy-yli
1
Fork 0
Code Issues 1 Pull Requests Releases Activity

Make mpmath optional dependency

Browse Source
This commit is contained in:
RunasSudo 2022-10-17 20:35:04 +11:00
parent 7a93355dab
commit ca521ae529
Signed by: RunasSudo
GPG Key ID: 7234E476BF21C61A
1 changed files with 42 additions and 31 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

73
yli/distributions.py
View File

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