Coverage for /usr/lib/python3/dist-packages/scipy/stats/_levy_stable/__init__.py: 16%

1# 

2 

3import warnings 

4from functools import partial 

5 

6import numpy as np 

7 

8from scipy import optimize 

9from scipy import integrate 

10from scipy.integrate._quadrature import _builtincoeffs 

11from scipy import interpolate 

12from scipy.interpolate import RectBivariateSpline 

13import scipy.special as sc 

14from scipy._lib._util import _lazywhere 

15from .._distn_infrastructure import rv_continuous, _ShapeInfo 

16from .._continuous_distns import uniform, expon, _norm_pdf, _norm_cdf 

17from .levyst import Nolan 

18from scipy._lib.doccer import inherit_docstring_from 

19 

20 

21__all__ = ["levy_stable", "levy_stable_gen", "pdf_from_cf_with_fft"] 

22 

23# Stable distributions are known for various parameterisations 

24# some being advantageous for numerical considerations and others 

25# useful due to their location/scale awareness. 

26# 

27# Here we follow [NO] convention (see the references in the docstring 

28# for levy_stable_gen below). 

29# 

30# S0 / Z0 / x0 (aka Zoleterav's M) 

31# S1 / Z1 / x1 

32# 

33# Where S* denotes parameterisation, Z* denotes standardized 

34# version where gamma = 1, delta = 0 and x* denotes variable. 

35# 

36# Scipy's original Stable was a random variate generator. It 

37# uses S1 and unfortunately is not a location/scale aware. 

38 

39 

40# default numerical integration tolerance 

41# used for epsrel in piecewise and both epsrel and epsabs in dni 

42# (epsabs needed in dni since weighted quad requires epsabs > 0) 

43_QUAD_EPS = 1.2e-14 

44 

45 

46def _Phi_Z0(alpha, t): 

47 return ( 

48 -np.tan(np.pi * alpha / 2) * (np.abs(t) ** (1 - alpha) - 1) 

49 if alpha != 1 

50 else -2.0 * np.log(np.abs(t)) / np.pi 

51 ) 

52 

53 

54def _Phi_Z1(alpha, t): 

55 return ( 

56 np.tan(np.pi * alpha / 2) 

57 if alpha != 1 

58 else -2.0 * np.log(np.abs(t)) / np.pi 

59 ) 

60 

61 

62def _cf(Phi, t, alpha, beta): 

63 """Characteristic function.""" 

64 return np.exp( 

65 -(np.abs(t) ** alpha) * (1 - 1j * beta * np.sign(t) * Phi(alpha, t)) 

66 ) 

67 

68 

69_cf_Z0 = partial(_cf, _Phi_Z0) 

70_cf_Z1 = partial(_cf, _Phi_Z1) 

71 

72 

73def _pdf_single_value_cf_integrate(Phi, x, alpha, beta, **kwds): 

74 """To improve DNI accuracy convert characteristic function in to real 

75 valued integral using Euler's formula, then exploit cosine symmetry to 

76 change limits to [0, inf). Finally use cosine addition formula to split 

77 into two parts that can be handled by weighted quad pack. 

78 """ 

79 quad_eps = kwds.get("quad_eps", _QUAD_EPS) 

80 

81 def integrand1(t): 

82 if t == 0: 

83 return 0 

84 return np.exp(-(t ** alpha)) * ( 

85 np.cos(beta * (t ** alpha) * Phi(alpha, t)) 

86 ) 

87 

88 def integrand2(t): 

89 if t == 0: 

90 return 0 

91 return np.exp(-(t ** alpha)) * ( 

92 np.sin(beta * (t ** alpha) * Phi(alpha, t)) 

93 ) 

94 

95 with np.errstate(invalid="ignore"): 

96 int1, *ret1 = integrate.quad( 

97 integrand1, 

98 0, 

99 np.inf, 

100 weight="cos", 

101 wvar=x, 

102 limit=1000, 

103 epsabs=quad_eps, 

104 epsrel=quad_eps, 

105 full_output=1, 

106 ) 

107 

108 int2, *ret2 = integrate.quad( 

109 integrand2, 

110 0, 

111 np.inf, 

112 weight="sin", 

113 wvar=x, 

114 limit=1000, 

115 epsabs=quad_eps, 

116 epsrel=quad_eps, 

117 full_output=1, 

118 ) 

119 

120 return (int1 + int2) / np.pi 

121 

122 

123_pdf_single_value_cf_integrate_Z0 = partial( 

124 _pdf_single_value_cf_integrate, _Phi_Z0 

125) 

126_pdf_single_value_cf_integrate_Z1 = partial( 

127 _pdf_single_value_cf_integrate, _Phi_Z1 

128) 

129 

130 

131def _nolan_round_difficult_input( 

132 x0, alpha, beta, zeta, x_tol_near_zeta, alpha_tol_near_one 

133): 

134 """Round difficult input values for Nolan's method in [NO].""" 

135 

136 # following Nolan's STABLE, 

137 # "1. When 0 < |alpha-1| < 0.005, the program has numerical problems 

138 # evaluating the pdf and cdf. The current version of the program sets 

139 # alpha=1 in these cases. This approximation is not bad in the S0 

140 # parameterization." 

141 if np.abs(alpha - 1) < alpha_tol_near_one: 

142 alpha = 1.0 

143 

144 # "2. When alpha=1 and |beta| < 0.005, the program has numerical 

145 # problems. The current version sets beta=0." 

146 # We seem to have addressed this through re-expression of g(theta) here 

147 

148 # "8. When |x0-beta*tan(pi*alpha/2)| is small, the 

149 # computations of the density and cumulative have numerical problems. 

150 # The program works around this by setting 

151 # z = beta*tan(pi*alpha/2) when 

152 # |z-beta*tan(pi*alpha/2)| < tol(5)*alpha**(1/alpha). 

153 # (The bound on the right is ad hoc, to get reasonable behavior 

154 # when alpha is small)." 

155 # where tol(5) = 0.5e-2 by default. 

156 # 

157 # We seem to have partially addressed this through re-expression of 

158 # g(theta) here, but it still needs to be used in some extreme cases. 

159 # Perhaps tol(5) = 0.5e-2 could be reduced for our implementation. 

160 if np.abs(x0 - zeta) < x_tol_near_zeta * alpha ** (1 / alpha): 

161 x0 = zeta 

162 

163 return x0, alpha, beta 

164 

165 

166def _pdf_single_value_piecewise_Z1(x, alpha, beta, **kwds): 

167 # convert from Nolan's S_1 (aka S) to S_0 (aka Zolaterev M) 

168 # parameterization 

169 

170 zeta = -beta * np.tan(np.pi * alpha / 2.0) 

171 x0 = x + zeta if alpha != 1 else x 

172 

173 return _pdf_single_value_piecewise_Z0(x0, alpha, beta, **kwds) 

174 

175 

176def _pdf_single_value_piecewise_Z0(x0, alpha, beta, **kwds): 

177 

178 quad_eps = kwds.get("quad_eps", _QUAD_EPS) 

179 x_tol_near_zeta = kwds.get("piecewise_x_tol_near_zeta", 0.005) 

180 alpha_tol_near_one = kwds.get("piecewise_alpha_tol_near_one", 0.005) 

181 

182 zeta = -beta * np.tan(np.pi * alpha / 2.0) 

183 x0, alpha, beta = _nolan_round_difficult_input( 

184 x0, alpha, beta, zeta, x_tol_near_zeta, alpha_tol_near_one 

185 ) 

186 

187 # some other known distribution pdfs / analytical cases 

188 # TODO: add more where possible with test coverage, 

189 # eg https://en.wikipedia.org/wiki/Stable_distribution#Other_analytic_cases 

190 if alpha == 2.0: 

191 # normal 

192 return _norm_pdf(x0 / np.sqrt(2)) / np.sqrt(2) 

193 elif alpha == 0.5 and beta == 1.0: 

194 # levy 

195 # since S(1/2, 1, gamma, delta; <x>) == 

196 # S(1/2, 1, gamma, gamma + delta; <x0>). 

197 _x = x0 + 1 

198 if _x <= 0: 

199 return 0 

200 

201 return 1 / np.sqrt(2 * np.pi * _x) / _x * np.exp(-1 / (2 * _x)) 

202 elif alpha == 0.5 and beta == 0.0 and x0 != 0: 

203 # analytical solution [HO] 

204 S, C = sc.fresnel([1 / np.sqrt(2 * np.pi * np.abs(x0))]) 

205 arg = 1 / (4 * np.abs(x0)) 

206 return ( 

207 np.sin(arg) * (0.5 - S[0]) + np.cos(arg) * (0.5 - C[0]) 

208 ) / np.sqrt(2 * np.pi * np.abs(x0) ** 3) 

209 elif alpha == 1.0 and beta == 0.0: 

210 # cauchy 

211 return 1 / (1 + x0 ** 2) / np.pi 

212 

213 return _pdf_single_value_piecewise_post_rounding_Z0( 

214 x0, alpha, beta, quad_eps 

215 ) 

216 

217 

218def _pdf_single_value_piecewise_post_rounding_Z0(x0, alpha, beta, quad_eps): 

219 """Calculate pdf using Nolan's methods as detailed in [NO]. 

220 """ 

221 

222 _nolan = Nolan(alpha, beta, x0) 

223 zeta = _nolan.zeta 

224 xi = _nolan.xi 

225 c2 = _nolan.c2 

226 g = _nolan.g 

227 

228 # handle Nolan's initial case logic 

229 if x0 == zeta: 

230 return ( 

231 sc.gamma(1 + 1 / alpha) 

232 * np.cos(xi) 

233 / np.pi 

234 / ((1 + zeta ** 2) ** (1 / alpha / 2)) 

235 ) 

236 elif x0 < zeta: 

237 return _pdf_single_value_piecewise_post_rounding_Z0( 

238 -x0, alpha, -beta, quad_eps 

239 ) 

240 

241 # following Nolan, we may now assume 

242 # x0 > zeta when alpha != 1 

243 # beta != 0 when alpha == 1 

244 

245 # spare calculating integral on null set 

246 # use isclose as macos has fp differences 

247 if np.isclose(-xi, np.pi / 2, rtol=1e-014, atol=1e-014): 

248 return 0.0 

249 

250 def integrand(theta): 

251 # limit any numerical issues leading to g_1 < 0 near theta limits 

252 g_1 = g(theta) 

253 if not np.isfinite(g_1) or g_1 < 0: 

254 g_1 = 0 

255 return g_1 * np.exp(-g_1) 

256 

257 with np.errstate(all="ignore"): 

258 peak = optimize.bisect( 

259 lambda t: g(t) - 1, -xi, np.pi / 2, xtol=quad_eps 

260 ) 

261 

262 # this integrand can be very peaked, so we need to force 

263 # QUADPACK to evaluate the function inside its support 

264 # 

265 

266 # lastly, we add additional samples at 

267 # ~exp(-100), ~exp(-10), ~exp(-5), ~exp(-1) 

268 # to improve QUADPACK's detection of rapidly descending tail behavior 

269 # (this choice is fairly ad hoc) 

270 tail_points = [ 

271 optimize.bisect(lambda t: g(t) - exp_height, -xi, np.pi / 2) 

272 for exp_height in [100, 10, 5] 

273 # exp_height = 1 is handled by peak 

274 ] 

275 intg_points = [0, peak] + tail_points 

276 intg, *ret = integrate.quad( 

277 integrand, 

278 -xi, 

279 np.pi / 2, 

280 points=intg_points, 

281 limit=100, 

282 epsrel=quad_eps, 

283 epsabs=0, 

284 full_output=1, 

285 ) 

286 

287 return c2 * intg 

288 

289 

290def _cdf_single_value_piecewise_Z1(x, alpha, beta, **kwds): 

291 # convert from Nolan's S_1 (aka S) to S_0 (aka Zolaterev M) 

292 # parameterization 

293 

294 zeta = -beta * np.tan(np.pi * alpha / 2.0) 

295 x0 = x + zeta if alpha != 1 else x 

296 

297 return _cdf_single_value_piecewise_Z0(x0, alpha, beta, **kwds) 

298 

299 

300def _cdf_single_value_piecewise_Z0(x0, alpha, beta, **kwds): 

301 

302 quad_eps = kwds.get("quad_eps", _QUAD_EPS) 

303 x_tol_near_zeta = kwds.get("piecewise_x_tol_near_zeta", 0.005) 

304 alpha_tol_near_one = kwds.get("piecewise_alpha_tol_near_one", 0.005) 

305 

306 zeta = -beta * np.tan(np.pi * alpha / 2.0) 

307 x0, alpha, beta = _nolan_round_difficult_input( 

308 x0, alpha, beta, zeta, x_tol_near_zeta, alpha_tol_near_one 

309 ) 

310 

311 # some other known distribution cdfs / analytical cases 

312 # TODO: add more where possible with test coverage, 

313 # eg https://en.wikipedia.org/wiki/Stable_distribution#Other_analytic_cases 

314 if alpha == 2.0: 

315 # normal 

316 return _norm_cdf(x0 / np.sqrt(2)) 

317 elif alpha == 0.5 and beta == 1.0: 

318 # levy 

319 # since S(1/2, 1, gamma, delta; <x>) == 

320 # S(1/2, 1, gamma, gamma + delta; <x0>). 

321 _x = x0 + 1 

322 if _x <= 0: 

323 return 0 

324 

325 return sc.erfc(np.sqrt(0.5 / _x)) 

326 elif alpha == 1.0 and beta == 0.0: 

327 # cauchy 

328 return 0.5 + np.arctan(x0) / np.pi 

329 

330 return _cdf_single_value_piecewise_post_rounding_Z0( 

331 x0, alpha, beta, quad_eps 

332 ) 

333 

334 

335def _cdf_single_value_piecewise_post_rounding_Z0(x0, alpha, beta, quad_eps): 

336 """Calculate cdf using Nolan's methods as detailed in [NO]. 

337 """ 

338 _nolan = Nolan(alpha, beta, x0) 

339 zeta = _nolan.zeta 

340 xi = _nolan.xi 

341 c1 = _nolan.c1 

342 # c2 = _nolan.c2 

343 c3 = _nolan.c3 

344 g = _nolan.g 

345 

346 # handle Nolan's initial case logic 

347 if (alpha == 1 and beta < 0) or x0 < zeta: 

348 # NOTE: Nolan's paper has a typo here! 

349 # He states F(x) = 1 - F(x, alpha, -beta), but this is clearly 

350 # incorrect since F(-infty) would be 1.0 in this case 

351 # Indeed, the alpha != 1, x0 < zeta case is correct here. 

352 return 1 - _cdf_single_value_piecewise_post_rounding_Z0( 

353 -x0, alpha, -beta, quad_eps 

354 ) 

355 elif x0 == zeta: 

356 return 0.5 - xi / np.pi 

357 

358 # following Nolan, we may now assume 

359 # x0 > zeta when alpha != 1 

360 # beta > 0 when alpha == 1 

361 

362 # spare calculating integral on null set 

363 # use isclose as macos has fp differences 

364 if np.isclose(-xi, np.pi / 2, rtol=1e-014, atol=1e-014): 

365 return c1 

366 

367 def integrand(theta): 

368 g_1 = g(theta) 

369 return np.exp(-g_1) 

370 

371 with np.errstate(all="ignore"): 

372 # shrink supports where required 

373 left_support = -xi 

374 right_support = np.pi / 2 

375 if alpha > 1: 

376 # integrand(t) monotonic 0 to 1 

377 if integrand(-xi) != 0.0: 

378 res = optimize.minimize( 

379 integrand, 

380 (-xi,), 

381 method="L-BFGS-B", 

382 bounds=[(-xi, np.pi / 2)], 

383 ) 

384 left_support = res.x[0] 

385 else: 

386 # integrand(t) monotonic 1 to 0 

387 if integrand(np.pi / 2) != 0.0: 

388 res = optimize.minimize( 

389 integrand, 

390 (np.pi / 2,), 

391 method="L-BFGS-B", 

392 bounds=[(-xi, np.pi / 2)], 

393 ) 

394 right_support = res.x[0] 

395 

396 intg, *ret = integrate.quad( 

397 integrand, 

398 left_support, 

399 right_support, 

400 points=[left_support, right_support], 

401 limit=100, 

402 epsrel=quad_eps, 

403 epsabs=0, 

404 full_output=1, 

405 ) 

406 

407 return c1 + c3 * intg 

408 

409 

410def _rvs_Z1(alpha, beta, size=None, random_state=None): 

411 """Simulate random variables using Nolan's methods as detailed in [NO]. 

412 """ 

413 

414 def alpha1func(alpha, beta, TH, aTH, bTH, cosTH, tanTH, W): 

415 return ( 

416 2 

417 / np.pi 

418 * ( 

419 (np.pi / 2 + bTH) * tanTH 

420 - beta * np.log((np.pi / 2 * W * cosTH) / (np.pi / 2 + bTH)) 

421 ) 

422 ) 

423 

424 def beta0func(alpha, beta, TH, aTH, bTH, cosTH, tanTH, W): 

425 return ( 

426 W 

427 / (cosTH / np.tan(aTH) + np.sin(TH)) 

428 * ((np.cos(aTH) + np.sin(aTH) * tanTH) / W) ** (1.0 / alpha) 

429 ) 

430 

431 def otherwise(alpha, beta, TH, aTH, bTH, cosTH, tanTH, W): 

432 # alpha is not 1 and beta is not 0 

433 val0 = beta * np.tan(np.pi * alpha / 2) 

434 th0 = np.arctan(val0) / alpha 

435 val3 = W / (cosTH / np.tan(alpha * (th0 + TH)) + np.sin(TH)) 

436 res3 = val3 * ( 

437 ( 

438 np.cos(aTH) 

439 + np.sin(aTH) * tanTH 

440 - val0 * (np.sin(aTH) - np.cos(aTH) * tanTH) 

441 ) 

442 / W 

443 ) ** (1.0 / alpha) 

444 return res3 

445 

446 def alphanot1func(alpha, beta, TH, aTH, bTH, cosTH, tanTH, W): 

447 res = _lazywhere( 

448 beta == 0, 

449 (alpha, beta, TH, aTH, bTH, cosTH, tanTH, W), 

450 beta0func, 

451 f2=otherwise, 

452 ) 

453 return res 

454 

455 alpha = np.broadcast_to(alpha, size) 

456 beta = np.broadcast_to(beta, size) 

457 TH = uniform.rvs( 

458 loc=-np.pi / 2.0, scale=np.pi, size=size, random_state=random_state 

459 ) 

460 W = expon.rvs(size=size, random_state=random_state) 

461 aTH = alpha * TH 

462 bTH = beta * TH 

463 cosTH = np.cos(TH) 

464 tanTH = np.tan(TH) 

465 res = _lazywhere( 

466 alpha == 1, 

467 (alpha, beta, TH, aTH, bTH, cosTH, tanTH, W), 

468 alpha1func, 

469 f2=alphanot1func, 

470 ) 

471 return res 

472 

473 

474def _fitstart_S0(data): 

475 alpha, beta, delta1, gamma = _fitstart_S1(data) 

476 

477 # Formulas for mapping parameters in S1 parameterization to 

478 # those in S0 parameterization can be found in [NO]. Note that 

479 # only delta changes. 

480 if alpha != 1: 

481 delta0 = delta1 + beta * gamma * np.tan(np.pi * alpha / 2.0) 

482 else: 

483 delta0 = delta1 + 2 * beta * gamma * np.log(gamma) / np.pi 

484 

485 return alpha, beta, delta0, gamma 

486 

487 

488def _fitstart_S1(data): 

489 # We follow McCullock 1986 method - Simple Consistent Estimators 

490 # of Stable Distribution Parameters 

491 

492 # fmt: off 

493 # Table III and IV 

494 nu_alpha_range = [2.439, 2.5, 2.6, 2.7, 2.8, 3, 3.2, 3.5, 4, 

495 5, 6, 8, 10, 15, 25] 

496 nu_beta_range = [0, 0.1, 0.2, 0.3, 0.5, 0.7, 1] 

497 

498 # table III - alpha = psi_1(nu_alpha, nu_beta) 

499 alpha_table = np.array([ 

500 [2.000, 2.000, 2.000, 2.000, 2.000, 2.000, 2.000], 

501 [1.916, 1.924, 1.924, 1.924, 1.924, 1.924, 1.924], 

502 [1.808, 1.813, 1.829, 1.829, 1.829, 1.829, 1.829], 

503 [1.729, 1.730, 1.737, 1.745, 1.745, 1.745, 1.745], 

504 [1.664, 1.663, 1.663, 1.668, 1.676, 1.676, 1.676], 

505 [1.563, 1.560, 1.553, 1.548, 1.547, 1.547, 1.547], 

506 [1.484, 1.480, 1.471, 1.460, 1.448, 1.438, 1.438], 

507 [1.391, 1.386, 1.378, 1.364, 1.337, 1.318, 1.318], 

508 [1.279, 1.273, 1.266, 1.250, 1.210, 1.184, 1.150], 

509 [1.128, 1.121, 1.114, 1.101, 1.067, 1.027, 0.973], 

510 [1.029, 1.021, 1.014, 1.004, 0.974, 0.935, 0.874], 

511 [0.896, 0.892, 0.884, 0.883, 0.855, 0.823, 0.769], 

512 [0.818, 0.812, 0.806, 0.801, 0.780, 0.756, 0.691], 

513 [0.698, 0.695, 0.692, 0.689, 0.676, 0.656, 0.597], 

514 [0.593, 0.590, 0.588, 0.586, 0.579, 0.563, 0.513]]).T 

515 # transpose because interpolation with `RectBivariateSpline` is with 

516 # `nu_beta` as `x` and `nu_alpha` as `y` 

517 

518 # table IV - beta = psi_2(nu_alpha, nu_beta) 

519 beta_table = np.array([ 

520 [0, 2.160, 1.000, 1.000, 1.000, 1.000, 1.000], 

521 [0, 1.592, 3.390, 1.000, 1.000, 1.000, 1.000], 

522 [0, 0.759, 1.800, 1.000, 1.000, 1.000, 1.000], 

523 [0, 0.482, 1.048, 1.694, 1.000, 1.000, 1.000], 

524 [0, 0.360, 0.760, 1.232, 2.229, 1.000, 1.000], 

525 [0, 0.253, 0.518, 0.823, 1.575, 1.000, 1.000], 

526 [0, 0.203, 0.410, 0.632, 1.244, 1.906, 1.000], 

527 [0, 0.165, 0.332, 0.499, 0.943, 1.560, 1.000], 

528 [0, 0.136, 0.271, 0.404, 0.689, 1.230, 2.195], 

529 [0, 0.109, 0.216, 0.323, 0.539, 0.827, 1.917], 

530 [0, 0.096, 0.190, 0.284, 0.472, 0.693, 1.759], 

531 [0, 0.082, 0.163, 0.243, 0.412, 0.601, 1.596], 

532 [0, 0.074, 0.147, 0.220, 0.377, 0.546, 1.482], 

533 [0, 0.064, 0.128, 0.191, 0.330, 0.478, 1.362], 

534 [0, 0.056, 0.112, 0.167, 0.285, 0.428, 1.274]]).T 

535 

536 # Table V and VII 

537 # These are ordered with decreasing `alpha_range`; so we will need to 

538 # reverse them as required by RectBivariateSpline. 

539 alpha_range = [2, 1.9, 1.8, 1.7, 1.6, 1.5, 1.4, 1.3, 1.2, 1.1, 

540 1, 0.9, 0.8, 0.7, 0.6, 0.5][::-1] 

541 beta_range = [0, 0.25, 0.5, 0.75, 1] 

542 

543 # Table V - nu_c = psi_3(alpha, beta) 

544 nu_c_table = np.array([ 

545 [1.908, 1.908, 1.908, 1.908, 1.908], 

546 [1.914, 1.915, 1.916, 1.918, 1.921], 

547 [1.921, 1.922, 1.927, 1.936, 1.947], 

548 [1.927, 1.930, 1.943, 1.961, 1.987], 

549 [1.933, 1.940, 1.962, 1.997, 2.043], 

550 [1.939, 1.952, 1.988, 2.045, 2.116], 

551 [1.946, 1.967, 2.022, 2.106, 2.211], 

552 [1.955, 1.984, 2.067, 2.188, 2.333], 

553 [1.965, 2.007, 2.125, 2.294, 2.491], 

554 [1.980, 2.040, 2.205, 2.435, 2.696], 

555 [2.000, 2.085, 2.311, 2.624, 2.973], 

556 [2.040, 2.149, 2.461, 2.886, 3.356], 

557 [2.098, 2.244, 2.676, 3.265, 3.912], 

558 [2.189, 2.392, 3.004, 3.844, 4.775], 

559 [2.337, 2.634, 3.542, 4.808, 6.247], 

560 [2.588, 3.073, 4.534, 6.636, 9.144]])[::-1].T 

561 # transpose because interpolation with `RectBivariateSpline` is with 

562 # `beta` as `x` and `alpha` as `y` 

563 

564 # Table VII - nu_zeta = psi_5(alpha, beta) 

565 nu_zeta_table = np.array([ 

566 [0, 0.000, 0.000, 0.000, 0.000], 

567 [0, -0.017, -0.032, -0.049, -0.064], 

568 [0, -0.030, -0.061, -0.092, -0.123], 

569 [0, -0.043, -0.088, -0.132, -0.179], 

570 [0, -0.056, -0.111, -0.170, -0.232], 

571 [0, -0.066, -0.134, -0.206, -0.283], 

572 [0, -0.075, -0.154, -0.241, -0.335], 

573 [0, -0.084, -0.173, -0.276, -0.390], 

574 [0, -0.090, -0.192, -0.310, -0.447], 

575 [0, -0.095, -0.208, -0.346, -0.508], 

576 [0, -0.098, -0.223, -0.380, -0.576], 

577 [0, -0.099, -0.237, -0.424, -0.652], 

578 [0, -0.096, -0.250, -0.469, -0.742], 

579 [0, -0.089, -0.262, -0.520, -0.853], 

580 [0, -0.078, -0.272, -0.581, -0.997], 

581 [0, -0.061, -0.279, -0.659, -1.198]])[::-1].T 

582 # fmt: on 

583 

584 psi_1 = RectBivariateSpline(nu_beta_range, nu_alpha_range, 

585 alpha_table, kx=1, ky=1, s=0) 

586 

587 def psi_1_1(nu_beta, nu_alpha): 

588 return psi_1(nu_beta, nu_alpha) \ 

589 if nu_beta > 0 else psi_1(-nu_beta, nu_alpha) 

590 

591 psi_2 = RectBivariateSpline(nu_beta_range, nu_alpha_range, 

592 beta_table, kx=1, ky=1, s=0) 

593 

594 def psi_2_1(nu_beta, nu_alpha): 

595 return psi_2(nu_beta, nu_alpha) \ 

596 if nu_beta > 0 else -psi_2(-nu_beta, nu_alpha) 

597 

598 phi_3 = RectBivariateSpline(beta_range, alpha_range, nu_c_table, 

599 kx=1, ky=1, s=0) 

600 

601 def phi_3_1(beta, alpha): 

602 return phi_3(beta, alpha) if beta > 0 else phi_3(-beta, alpha) 

603 

604 phi_5 = RectBivariateSpline(beta_range, alpha_range, nu_zeta_table, 

605 kx=1, ky=1, s=0) 

606 

607 def phi_5_1(beta, alpha): 

608 return phi_5(beta, alpha) if beta > 0 else -phi_5(-beta, alpha) 

609 

610 # quantiles 

611 p05 = np.percentile(data, 5) 

612 p50 = np.percentile(data, 50) 

613 p95 = np.percentile(data, 95) 

614 p25 = np.percentile(data, 25) 

615 p75 = np.percentile(data, 75) 

616 

617 nu_alpha = (p95 - p05) / (p75 - p25) 

618 nu_beta = (p95 + p05 - 2 * p50) / (p95 - p05) 

619 

620 if nu_alpha >= 2.439: 

621 eps = np.finfo(float).eps 

622 alpha = np.clip(psi_1_1(nu_beta, nu_alpha)[0, 0], eps, 2.) 

623 beta = np.clip(psi_2_1(nu_beta, nu_alpha)[0, 0], -1.0, 1.0) 

624 else: 

625 alpha = 2.0 

626 beta = np.sign(nu_beta) 

627 c = (p75 - p25) / phi_3_1(beta, alpha)[0, 0] 

628 zeta = p50 + c * phi_5_1(beta, alpha)[0, 0] 

629 delta = zeta-beta*c*np.tan(np.pi*alpha/2.) if alpha != 1. else zeta 

630 

631 return (alpha, beta, delta, c) 

632 

633 

634class levy_stable_gen(rv_continuous): 

635 r"""A Levy-stable continuous random variable. 

636 

637 %(before_notes)s 

638 

639 See Also 

640 -------- 

641 levy, levy_l, cauchy, norm 

642 

643 Notes 

644 ----- 

645 The distribution for `levy_stable` has characteristic function: 

646 

647 .. math:: 

648 

649 \varphi(t, \alpha, \beta, c, \mu) = 

650 e^{it\mu -|ct|^{\alpha}(1-i\beta\operatorname{sign}(t)\Phi(\alpha, t))} 

651 

652 where two different parameterizations are supported. The first :math:`S_1`: 

653 

654 .. math:: 

655 

656 \Phi = \begin{cases} 

657 \tan \left({\frac {\pi \alpha }{2}}\right)&\alpha \neq 1\\ 

658 -{\frac {2}{\pi }}\log |t|&\alpha =1 

659 \end{cases} 

660 

661 The second :math:`S_0`: 

662 

663 .. math:: 

664 

665 \Phi = \begin{cases} 

666 -\tan \left({\frac {\pi \alpha }{2}}\right)(|ct|^{1-\alpha}-1) 

667 &\alpha \neq 1\\ 

668 -{\frac {2}{\pi }}\log |ct|&\alpha =1 

669 \end{cases} 

670 

671 

672 The probability density function for `levy_stable` is: 

673 

674 .. math:: 

675 

676 f(x) = \frac{1}{2\pi}\int_{-\infty}^\infty \varphi(t)e^{-ixt}\,dt 

677 

678 where :math:`-\infty < t < \infty`. This integral does not have a known 

679 closed form. 

680 

681 `levy_stable` generalizes several distributions. Where possible, they 

682 should be used instead. Specifically, when the shape parameters 

683 assume the values in the table below, the corresponding equivalent 

684 distribution should be used. 

685 

686 ========= ======== =========== 

687 ``alpha`` ``beta`` Equivalent 

688 ========= ======== =========== 

689 1/2 -1 `levy_l` 

690 1/2 1 `levy` 

691 1 0 `cauchy` 

692 2 any `norm` (with ``scale=sqrt(2)``) 

693 ========= ======== =========== 

694 

695 Evaluation of the pdf uses Nolan's piecewise integration approach with the 

696 Zolotarev :math:`M` parameterization by default. There is also the option 

697 to use direct numerical integration of the standard parameterization of the 

698 characteristic function or to evaluate by taking the FFT of the 

699 characteristic function. 

700 

701 The default method can changed by setting the class variable 

702 ``levy_stable.pdf_default_method`` to one of 'piecewise' for Nolan's 

703 approach, 'dni' for direct numerical integration, or 'fft-simpson' for the 

704 FFT based approach. For the sake of backwards compatibility, the methods 

705 'best' and 'zolotarev' are equivalent to 'piecewise' and the method 

706 'quadrature' is equivalent to 'dni'. 

707 

708 The parameterization can be changed by setting the class variable 

709 ``levy_stable.parameterization`` to either 'S0' or 'S1'. 

710 The default is 'S1'. 

711 

712 To improve performance of piecewise and direct numerical integration one 

713 can specify ``levy_stable.quad_eps`` (defaults to 1.2e-14). This is used 

714 as both the absolute and relative quadrature tolerance for direct numerical 

715 integration and as the relative quadrature tolerance for the piecewise 

716 method. One can also specify ``levy_stable.piecewise_x_tol_near_zeta`` 

717 (defaults to 0.005) for how close x is to zeta before it is considered the 

718 same as x [NO]. The exact check is 

719 ``abs(x0 - zeta) < piecewise_x_tol_near_zeta*alpha**(1/alpha)``. One can 

720 also specify ``levy_stable.piecewise_alpha_tol_near_one`` (defaults to 

721 0.005) for how close alpha is to 1 before being considered equal to 1. 

722 

723 To increase accuracy of FFT calculation one can specify 

724 ``levy_stable.pdf_fft_grid_spacing`` (defaults to 0.001) and 

725 ``pdf_fft_n_points_two_power`` (defaults to None which means a value is 

726 calculated that sufficiently covers the input range). 

727 

728 Further control over FFT calculation is available by setting 

729 ``pdf_fft_interpolation_degree`` (defaults to 3) for spline order and 

730 ``pdf_fft_interpolation_level`` for determining the number of points to use 

731 in the Newton-Cotes formula when approximating the characteristic function 

732 (considered experimental). 

733 

734 Evaluation of the cdf uses Nolan's piecewise integration approach with the 

735 Zolatarev :math:`S_0` parameterization by default. There is also the option 

736 to evaluate through integration of an interpolated spline of the pdf 

737 calculated by means of the FFT method. The settings affecting FFT 

738 calculation are the same as for pdf calculation. The default cdf method can 

739 be changed by setting ``levy_stable.cdf_default_method`` to either 

740 'piecewise' or 'fft-simpson'. For cdf calculations the Zolatarev method is 

741 superior in accuracy, so FFT is disabled by default. 

742 

743 Fitting estimate uses quantile estimation method in [MC]. MLE estimation of 

744 parameters in fit method uses this quantile estimate initially. Note that 

745 MLE doesn't always converge if using FFT for pdf calculations; this will be 

746 the case if alpha <= 1 where the FFT approach doesn't give good 

747 approximations. 

748 

749 Any non-missing value for the attribute 

750 ``levy_stable.pdf_fft_min_points_threshold`` will set 

751 ``levy_stable.pdf_default_method`` to 'fft-simpson' if a valid 

752 default method is not otherwise set. 

753 

754 

755 

756 .. warning:: 

757 

758 For pdf calculations FFT calculation is considered experimental. 

759 

760 For cdf calculations FFT calculation is considered experimental. Use 

761 Zolatarev's method instead (default). 

762 

763 %(after_notes)s 

764 

765 References 

766 ---------- 

767 .. [MC] McCulloch, J., 1986. Simple consistent estimators of stable 

768 distribution parameters. Communications in Statistics - Simulation and 

769 Computation 15, 11091136. 

770 .. [WZ] Wang, Li and Zhang, Ji-Hong, 2008. Simpson's rule based FFT method 

771 to compute densities of stable distribution. 

772 .. [NO] Nolan, J., 1997. Numerical Calculation of Stable Densities and 

773 distributions Functions. 

774 .. [HO] Hopcraft, K. I., Jakeman, E., Tanner, R. M. J., 1999. Lévy random 

775 walks with fluctuating step number and multiscale behavior. 

776 

777 %(example)s 

778 

779 """ 

780 # Configurable options as class variables 

781 # (accesible from self by attribute lookup). 

782 parameterization = "S1" 

783 pdf_default_method = "piecewise" 

784 cdf_default_method = "piecewise" 

785 quad_eps = _QUAD_EPS 

786 piecewise_x_tol_near_zeta = 0.005 

787 piecewise_alpha_tol_near_one = 0.005