Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/stats/_levy_stable/__init__.py: 16%

1# -*- coding: utf-8 -*- 

2# 

3 

4import warnings 

5from functools import partial 

6 

7import numpy as np 

8 

9from scipy import optimize 

10from scipy import integrate 

11from scipy.integrate._quadrature import _builtincoeffs 

12from scipy import interpolate 

13from scipy.interpolate import RectBivariateSpline 

14import scipy.special as sc 

15from scipy._lib._util import _lazywhere 

16from .._distn_infrastructure import rv_continuous, _ShapeInfo 

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

18from .levyst import Nolan 

19from scipy._lib.doccer import inherit_docstring_from 

20 

21 

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

23 

24# Stable distributions are known for various parameterisations 

25# some being advantageous for numerical considerations and others 

26# useful due to their location/scale awareness. 

27# 

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

29# for levy_stable_gen below). 

30# 

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

32# S1 / Z1 / x1 

33# 

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

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

36# 

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

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

39 

40 

41# default numerical integration tolerance 

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

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

44_QUAD_EPS = 1.2e-14 

45 

46 

47def _Phi_Z0(alpha, t): 

48 return ( 

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

50 if alpha != 1 

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

52 ) 

53 

54 

55def _Phi_Z1(alpha, t): 

56 return ( 

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

58 if alpha != 1 

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

60 ) 

61 

62 

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

64 """Characteristic function.""" 

65 return np.exp( 

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

67 ) 

68 

69 

70_cf_Z0 = partial(_cf, _Phi_Z0) 

71_cf_Z1 = partial(_cf, _Phi_Z1) 

72 

73 

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

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

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

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

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

79 """ 

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

81 

82 def integrand1(t): 

83 if t == 0: 

84 return 0 

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

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

87 ) 

88 

89 def integrand2(t): 

90 if t == 0: 

91 return 0 

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

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

94 ) 

95 

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

97 int1, *ret1 = integrate.quad( 

98 integrand1, 

99 0, 

100 np.inf, 

101 weight="cos", 

102 wvar=x, 

103 limit=1000, 

104 epsabs=quad_eps, 

105 epsrel=quad_eps, 

106 full_output=1, 

107 ) 

108 

109 int2, *ret2 = integrate.quad( 

110 integrand2, 

111 0, 

112 np.inf, 

113 weight="sin", 

114 wvar=x, 

115 limit=1000, 

116 epsabs=quad_eps, 

117 epsrel=quad_eps, 

118 full_output=1, 

119 ) 

120 

121 return (int1 + int2) / np.pi 

122 

123 

124_pdf_single_value_cf_integrate_Z0 = partial( 

125 _pdf_single_value_cf_integrate, _Phi_Z0 

126) 

127_pdf_single_value_cf_integrate_Z1 = partial( 

128 _pdf_single_value_cf_integrate, _Phi_Z1 

129) 

130 

131 

132def _nolan_round_difficult_input( 

133 x0, alpha, beta, zeta, x_tol_near_zeta, alpha_tol_near_one 

134): 

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

136 

137 # following Nolan's STABLE, 

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

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

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

141 # parameterization." 

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

143 alpha = 1.0 

144 

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

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

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

148 

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

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

151 # The program works around this by setting 

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

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

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

155 # when alpha is small)." 

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

157 # 

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

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

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

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

162 x0 = zeta 

163 

164 return x0, alpha, beta 

165 

166 

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

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

169 # parameterization 

170 

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

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

173 

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

175 

176 

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

178 

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

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

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

182 

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

184 x0, alpha, beta = _nolan_round_difficult_input( 

185 x0, alpha, beta, zeta, x_tol_near_zeta, alpha_tol_near_one 

186 ) 

187 

188 # some other known distribution pdfs / analytical cases 

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

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

191 if alpha == 2.0: 

192 # normal 

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

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

195 # levy 

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

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

198 _x = x0 + 1 

199 if _x <= 0: 

200 return 0 

201 

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

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

204 # analytical solution [HO] 

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

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

207 return ( 

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

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

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

211 # cauchy 

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

213 

214 return _pdf_single_value_piecewise_post_rounding_Z0( 

215 x0, alpha, beta, quad_eps 

216 ) 

217 

218 

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

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

221 """ 

222 

223 _nolan = Nolan(alpha, beta, x0) 

224 zeta = _nolan.zeta 

225 xi = _nolan.xi 

226 c2 = _nolan.c2 

227 g = _nolan.g 

228 

229 # handle Nolan's initial case logic 

230 if x0 == zeta: 

231 return ( 

232 sc.gamma(1 + 1 / alpha) 

233 * np.cos(xi) 

234 / np.pi 

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

236 ) 

237 elif x0 < zeta: 

238 return _pdf_single_value_piecewise_post_rounding_Z0( 

239 -x0, alpha, -beta, quad_eps 

240 ) 

241 

242 # following Nolan, we may now assume 

243 # x0 > zeta when alpha != 1 

244 # beta != 0 when alpha == 1 

245 

246 # spare calculating integral on null set 

247 # use isclose as macos has fp differences 

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

249 return 0.0 

250 

251 def integrand(theta): 

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

253 g_1 = g(theta) 

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

255 g_1 = 0 

256 return g_1 * np.exp(-g_1) 

257 

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

259 peak = optimize.bisect( 

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

261 ) 

262 

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

264 # QUADPACK to evaluate the function inside its support 

265 # 

266 

267 # lastly, we add additional samples at 

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

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

270 # (this choice is fairly ad hoc) 

271 tail_points = [ 

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

273 for exp_height in [100, 10, 5] 

274 # exp_height = 1 is handled by peak 

275 ] 

276 intg_points = [0, peak] + tail_points 

277 intg, *ret = integrate.quad( 

278 integrand, 

279 -xi, 

280 np.pi / 2, 

281 points=intg_points, 

282 limit=100, 

283 epsrel=quad_eps, 

284 epsabs=0, 

285 full_output=1, 

286 ) 

287 

288 return c2 * intg 

289 

290 

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

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

293 # parameterization 

294 

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

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

297 

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

299 

300 

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

302 

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

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

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

306 

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

308 x0, alpha, beta = _nolan_round_difficult_input( 

309 x0, alpha, beta, zeta, x_tol_near_zeta, alpha_tol_near_one 

310 ) 

311 

312 # some other known distribution cdfs / analytical cases 

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

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

315 if alpha == 2.0: 

316 # normal 

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

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

319 # levy 

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

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

322 _x = x0 + 1 

323 if _x <= 0: 

324 return 0 

325 

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

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

328 # cauchy 

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

330 

331 return _cdf_single_value_piecewise_post_rounding_Z0( 

332 x0, alpha, beta, quad_eps 

333 ) 

334 

335 

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

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

338 """ 

339 _nolan = Nolan(alpha, beta, x0) 

340 zeta = _nolan.zeta 

341 xi = _nolan.xi 

342 c1 = _nolan.c1 

343 # c2 = _nolan.c2 

344 c3 = _nolan.c3 

345 g = _nolan.g 

346 

347 # handle Nolan's initial case logic 

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

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

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

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

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

353 return 1 - _cdf_single_value_piecewise_post_rounding_Z0( 

354 -x0, alpha, -beta, quad_eps 

355 ) 

356 elif x0 == zeta: 

357 return 0.5 - xi / np.pi 

358 

359 # following Nolan, we may now assume 

360 # x0 > zeta when alpha != 1 

361 # beta > 0 when alpha == 1 

362 

363 # spare calculating integral on null set 

364 # use isclose as macos has fp differences 

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

366 return c1 

367 

368 def integrand(theta): 

369 g_1 = g(theta) 

370 return np.exp(-g_1) 

371 

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

373 # shrink supports where required 

374 left_support = -xi 

375 right_support = np.pi / 2 

376 if alpha > 1: 

377 # integrand(t) monotonic 0 to 1 

378 if integrand(-xi) != 0.0: 

379 res = optimize.minimize( 

380 integrand, 

381 (-xi,), 

382 method="L-BFGS-B", 

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

384 ) 

385 left_support = res.x[0] 

386 else: 

387 # integrand(t) monotonic 1 to 0 

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

389 res = optimize.minimize( 

390 integrand, 

391 (np.pi / 2,), 

392 method="L-BFGS-B", 

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

394 ) 

395 right_support = res.x[0] 

396 

397 intg, *ret = integrate.quad( 

398 integrand, 

399 left_support, 

400 right_support, 

401 points=[left_support, right_support], 

402 limit=100, 

403 epsrel=quad_eps, 

404 epsabs=0, 

405 full_output=1, 

406 ) 

407 

408 return c1 + c3 * intg 

409 

410 

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

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

413 """ 

414 

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

416 return ( 

417 2 

418 / np.pi 

419 * ( 

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

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

422 ) 

423 ) 

424 

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

426 return ( 

427 W 

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

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

430 ) 

431 

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

433 # alpha is not 1 and beta is not 0 

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

435 th0 = np.arctan(val0) / alpha 

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

437 res3 = val3 * ( 

438 ( 

439 np.cos(aTH) 

440 + np.sin(aTH) * tanTH 

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

442 ) 

443 / W 

444 ) ** (1.0 / alpha) 

445 return res3 

446 

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

448 res = _lazywhere( 

449 beta == 0, 

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

451 beta0func, 

452 f2=otherwise, 

453 ) 

454 return res 

455 

456 alpha = np.broadcast_to(alpha, size) 

457 beta = np.broadcast_to(beta, size) 

458 TH = uniform.rvs( 

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

460 ) 

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

462 aTH = alpha * TH 

463 bTH = beta * TH 

464 cosTH = np.cos(TH) 

465 tanTH = np.tan(TH) 

466 res = _lazywhere( 

467 alpha == 1, 

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

469 alpha1func, 

470 f2=alphanot1func, 

471 ) 

472 return res 

473 

474 

475def _fitstart_S0(data): 

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

477 

478 # Formulas for mapping parameters in S1 parameterization to 

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

480 # only delta changes. 

481 if alpha != 1: 

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

483 else: 

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

485 

486 return alpha, beta, delta0, gamma 

487 

488 

489def _fitstart_S1(data): 

490 # We follow McCullock 1986 method - Simple Consistent Estimators 

491 # of Stable Distribution Parameters 

492 

493 # fmt: off 

494 # Table III and IV 

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

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

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

498 

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

500 alpha_table = np.array([ 

501 [2.000, 2.000, 2.000, 2.000, 2.000, 2.000, 2.000], 

502 [1.916, 1.924, 1.924, 1.924, 1.924, 1.924, 1.924], 

503 [1.808, 1.813, 1.829, 1.829, 1.829, 1.829, 1.829], 

504 [1.729, 1.730, 1.737, 1.745, 1.745, 1.745, 1.745], 

505 [1.664, 1.663, 1.663, 1.668, 1.676, 1.676, 1.676], 

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

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

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

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

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

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

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

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

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

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

516 # transpose because interpolation with `RectBivariateSpline` is with 

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

518 

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

520 beta_table = np.array([ 

521 [0, 2.160, 1.000, 1.000, 1.000, 1.000, 1.000], 

522 [0, 1.592, 3.390, 1.000, 1.000, 1.000, 1.000], 

523 [0, 0.759, 1.800, 1.000, 1.000, 1.000, 1.000], 

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

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

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

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

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

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

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

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

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

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

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

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

536 

537 # Table V and VII 

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

539 # reverse them as required by RectBivariateSpline. 

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

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

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

543 

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

545 nu_c_table = np.array([ 

546 [1.908, 1.908, 1.908, 1.908, 1.908], 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

562 # transpose because interpolation with `RectBivariateSpline` is with 

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

564 

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

566 nu_zeta_table = np.array([ 

567 [0, 0.000, 0.000, 0.000, 0.000], 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

583 # fmt: on 

584 

585 psi_1 = RectBivariateSpline(nu_beta_range, nu_alpha_range, 

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

587 

588 def psi_1_1(nu_beta, nu_alpha): 

589 return psi_1(nu_beta, nu_alpha) \ 

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

591 

592 psi_2 = RectBivariateSpline(nu_beta_range, nu_alpha_range, 

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

594 

595 def psi_2_1(nu_beta, nu_alpha): 

596 return psi_2(nu_beta, nu_alpha) \ 

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

598 

599 phi_3 = RectBivariateSpline(beta_range, alpha_range, nu_c_table, 

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

601 

602 def phi_3_1(beta, alpha): 

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

604 

605 phi_5 = RectBivariateSpline(beta_range, alpha_range, nu_zeta_table, 

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

607 

608 def phi_5_1(beta, alpha): 

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

610 

611 # quantiles 

612 p05 = np.percentile(data, 5) 

613 p50 = np.percentile(data, 50) 

614 p95 = np.percentile(data, 95) 

615 p25 = np.percentile(data, 25) 

616 p75 = np.percentile(data, 75) 

617 

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

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

620 

621 if nu_alpha >= 2.439: 

622 eps = np.finfo(float).eps 

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

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

625 else: 

626 alpha = 2.0 

627 beta = np.sign(nu_beta) 

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

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

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

631 

632 return (alpha, beta, delta, c) 

633 

634 

635class levy_stable_gen(rv_continuous): 

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

637 

638 %(before_notes)s 

639 

640 See Also 

641 -------- 

642 levy, levy_l, cauchy, norm 

643 

644 Notes 

645 ----- 

646 The distribution for `levy_stable` has characteristic function: 

647 

648 .. math:: 

649 

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

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

652 

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

654 

655 .. math:: 

656 

657 \Phi = \begin{cases} 

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

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

660 \end{cases} 

661 

662 The second :math:`S_0`: 

663 

664 .. math:: 

665 

666 \Phi = \begin{cases} 

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

668 &\alpha \neq 1\\ 

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

670 \end{cases} 

671 

672 

673 The probability density function for `levy_stable` is: 

674 

675 .. math:: 

676 

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

678 

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

680 closed form. 

681 

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

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

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

685 distribution should be used. 

686 

687 ========= ======== =========== 

688 ``alpha`` ``beta`` Equivalent 

689 ========= ======== =========== 

690 1/2 -1 `levy_l` 

691 1/2 1 `levy` 

692 1 0 `cauchy` 

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

694 ========= ======== =========== 

695 

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

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

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

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

700 characteristic function. 

701 

702 The default method can changed by setting the class variable 

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

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

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

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

707 'quadrature' is equivalent to 'dni'. 

708 

709 The parameterization can be changed by setting the class variable 

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

711 The default is 'S1'. 

712 

713 To improve performance of piecewise and direct numerical integration one 

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

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

716 integration and as the relative quadrature tolerance for the piecewise 

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

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

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

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

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

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

723 

724 To increase accuracy of FFT calculation one can specify 

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

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

727 calculated that sufficiently covers the input range). 

728 

729 Further control over FFT calculation is available by setting 

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

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

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

733 (considered experimental). 

734 

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

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

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

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

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

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

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

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

743 

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

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

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

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

748 approximations. 

749 

750 Any non-missing value for the attribute 

751 ``levy_stable.pdf_fft_min_points_threshold`` will set 

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

753 default method is not otherwise set. 

754 

755 

756 

757 .. warning:: 

758 

759 For pdf calculations FFT calculation is considered experimental. 

760 

761 For cdf calculations FFT calculation is considered experimental. Use 

762 Zolatarev's method instead (default). 

763 

764 %(after_notes)s 

765 

766 References 

767 ---------- 

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

769 distribution parameters. Communications in Statistics - Simulation and 

770 Computation 15, 11091136. 

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

772 to compute densities of stable distribution. 

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

774 distributions Functions. 

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

776 walks with fluctuating step number and multiscale behavior. 

777 

778 %(example)s 

779 

780 """ 

781 # Configurable options as class variables 

782 # (accesible from self by attribute lookup). 

783 parameterization = "S1" 

784 pdf_default_method = "piecewise" 

785 cdf_default_method = "piecewise" 

786 quad_eps = _QUAD_EPS 

787 piecewise_x_tol_near_zeta = 0.005 

788 piecewise_alpha_tol_near_one = 0.005