Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/stats/_tukeylambda_stats.py: 25%
53 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1import numpy as np
2from numpy import poly1d
3from scipy.special import beta
6# The following code was used to generate the Pade coefficients for the
7# Tukey Lambda variance function. Version 0.17 of mpmath was used.
8#---------------------------------------------------------------------------
9# import mpmath as mp
10#
11# mp.mp.dps = 60
12#
13# one = mp.mpf(1)
14# two = mp.mpf(2)
15#
16# def mpvar(lam):
17# if lam == 0:
18# v = mp.pi**2 / three
19# else:
20# v = (two / lam**2) * (one / (one + two*lam) -
21# mp.beta(lam + one, lam + one))
22# return v
23#
24# t = mp.taylor(mpvar, 0, 8)
25# p, q = mp.pade(t, 4, 4)
26# print("p =", [mp.fp.mpf(c) for c in p])
27# print("q =", [mp.fp.mpf(c) for c in q])
28#---------------------------------------------------------------------------
30# Pade coefficients for the Tukey Lambda variance function.
31_tukeylambda_var_pc = [3.289868133696453, 0.7306125098871127,
32 -0.5370742306855439, 0.17292046290190008,
33 -0.02371146284628187]
34_tukeylambda_var_qc = [1.0, 3.683605511659861, 4.184152498888124,
35 1.7660926747377275, 0.2643989311168465]
37# numpy.poly1d instances for the numerator and denominator of the
38# Pade approximation to the Tukey Lambda variance.
39_tukeylambda_var_p = poly1d(_tukeylambda_var_pc[::-1])
40_tukeylambda_var_q = poly1d(_tukeylambda_var_qc[::-1])
43def tukeylambda_variance(lam):
44 """Variance of the Tukey Lambda distribution.
46 Parameters
47 ----------
48 lam : array_like
49 The lambda values at which to compute the variance.
51 Returns
52 -------
53 v : ndarray
54 The variance. For lam < -0.5, the variance is not defined, so
55 np.nan is returned. For lam = 0.5, np.inf is returned.
57 Notes
58 -----
59 In an interval around lambda=0, this function uses the [4,4] Pade
60 approximation to compute the variance. Otherwise it uses the standard
61 formula (https://en.wikipedia.org/wiki/Tukey_lambda_distribution). The
62 Pade approximation is used because the standard formula has a removable
63 discontinuity at lambda = 0, and does not produce accurate numerical
64 results near lambda = 0.
65 """
66 lam = np.asarray(lam)
67 shp = lam.shape
68 lam = np.atleast_1d(lam).astype(np.float64)
70 # For absolute values of lam less than threshold, use the Pade
71 # approximation.
72 threshold = 0.075
74 # Play games with masks to implement the conditional evaluation of
75 # the distribution.
76 # lambda < -0.5: var = nan
77 low_mask = lam < -0.5
78 # lambda == -0.5: var = inf
79 neghalf_mask = lam == -0.5
80 # abs(lambda) < threshold: use Pade approximation
81 small_mask = np.abs(lam) < threshold
82 # else the "regular" case: use the explicit formula.
83 reg_mask = ~(low_mask | neghalf_mask | small_mask)
85 # Get the 'lam' values for the cases where they are needed.
86 small = lam[small_mask]
87 reg = lam[reg_mask]
89 # Compute the function for each case.
90 v = np.empty_like(lam)
91 v[low_mask] = np.nan
92 v[neghalf_mask] = np.inf
93 if small.size > 0:
94 # Use the Pade approximation near lambda = 0.
95 v[small_mask] = _tukeylambda_var_p(small) / _tukeylambda_var_q(small)
96 if reg.size > 0:
97 v[reg_mask] = (2.0 / reg**2) * (1.0 / (1.0 + 2 * reg) -
98 beta(reg + 1, reg + 1))
99 v.shape = shp
100 return v
103# The following code was used to generate the Pade coefficients for the
104# Tukey Lambda kurtosis function. Version 0.17 of mpmath was used.
105#---------------------------------------------------------------------------
106# import mpmath as mp
107#
108# mp.mp.dps = 60
109#
110# one = mp.mpf(1)
111# two = mp.mpf(2)
112# three = mp.mpf(3)
113# four = mp.mpf(4)
114#
115# def mpkurt(lam):
116# if lam == 0:
117# k = mp.mpf(6)/5
118# else:
119# numer = (one/(four*lam+one) - four*mp.beta(three*lam+one, lam+one) +
120# three*mp.beta(two*lam+one, two*lam+one))
121# denom = two*(one/(two*lam+one) - mp.beta(lam+one,lam+one))**2
122# k = numer / denom - three
123# return k
124#
125# # There is a bug in mpmath 0.17: when we use the 'method' keyword of the
126# # taylor function and we request a degree 9 Taylor polynomial, we actually
127# # get degree 8.
128# t = mp.taylor(mpkurt, 0, 9, method='quad', radius=0.01)
129# t = [mp.chop(c, tol=1e-15) for c in t]
130# p, q = mp.pade(t, 4, 4)
131# print("p =", [mp.fp.mpf(c) for c in p])
132# print("q =", [mp.fp.mpf(c) for c in q])
133#---------------------------------------------------------------------------
135# Pade coefficients for the Tukey Lambda kurtosis function.
136_tukeylambda_kurt_pc = [1.2, -5.853465139719495, -22.653447381131077,
137 0.20601184383406815, 4.59796302262789]
138_tukeylambda_kurt_qc = [1.0, 7.171149192233599, 12.96663094361842,
139 0.43075235247853005, -2.789746758009912]
141# numpy.poly1d instances for the numerator and denominator of the
142# Pade approximation to the Tukey Lambda kurtosis.
143_tukeylambda_kurt_p = poly1d(_tukeylambda_kurt_pc[::-1])
144_tukeylambda_kurt_q = poly1d(_tukeylambda_kurt_qc[::-1])
147def tukeylambda_kurtosis(lam):
148 """Kurtosis of the Tukey Lambda distribution.
150 Parameters
151 ----------
152 lam : array_like
153 The lambda values at which to compute the variance.
155 Returns
156 -------
157 v : ndarray
158 The variance. For lam < -0.25, the variance is not defined, so
159 np.nan is returned. For lam = 0.25, np.inf is returned.
161 """
162 lam = np.asarray(lam)
163 shp = lam.shape
164 lam = np.atleast_1d(lam).astype(np.float64)
166 # For absolute values of lam less than threshold, use the Pade
167 # approximation.
168 threshold = 0.055
170 # Use masks to implement the conditional evaluation of the kurtosis.
171 # lambda < -0.25: kurtosis = nan
172 low_mask = lam < -0.25
173 # lambda == -0.25: kurtosis = inf
174 negqrtr_mask = lam == -0.25
175 # lambda near 0: use Pade approximation
176 small_mask = np.abs(lam) < threshold
177 # else the "regular" case: use the explicit formula.
178 reg_mask = ~(low_mask | negqrtr_mask | small_mask)
180 # Get the 'lam' values for the cases where they are needed.
181 small = lam[small_mask]
182 reg = lam[reg_mask]
184 # Compute the function for each case.
185 k = np.empty_like(lam)
186 k[low_mask] = np.nan
187 k[negqrtr_mask] = np.inf
188 if small.size > 0:
189 k[small_mask] = _tukeylambda_kurt_p(small) / _tukeylambda_kurt_q(small)
190 if reg.size > 0:
191 numer = (1.0 / (4 * reg + 1) - 4 * beta(3 * reg + 1, reg + 1) +
192 3 * beta(2 * reg + 1, 2 * reg + 1))
193 denom = 2 * (1.0/(2 * reg + 1) - beta(reg + 1, reg + 1))**2
194 k[reg_mask] = numer / denom - 3
196 # The return value will be a numpy array; resetting the shape ensures that
197 # if `lam` was a scalar, the return value is a 0-d array.
198 k.shape = shp
199 return k