DOC/ENH: Explain limitations of stats.levy_stable more explicitly

[jcalbert]

In response to #8263 - makes the documentation of stats.levy_stable more descriptive.

[ev-br]

Cross-ref #7374

[chrisb83]

Hi,
having these comments is definitely an improvement. I added a few remarks...

[chrisb83]

"function cannot be expressed analytically" --> "function cannot be expressed analytically (except in some special cases)"

[jcalbert]

I agree

[chrisb83]

#7374 hast been merged in the meantime. The table with the special cases would still be a useful addition. @jcalbert do you want to rebase and update your changes?

[jcalbert]

Sure, I'll do that early next week. About #7374 - is the new _pdf() sufficient for all purposes, or is it still preferable to use other distributions (e.g. cauchy for alpha=1, beta=0) when possible?

[chrisb83]

Thanks for the update. I don't know the implementation of 7374, my guess is that the special cases use simpler methods.

[jcalbert]

OK, I think it's good but my rebase skills are shaky so we'll see.

[chrisb83]

Once the rebase is successful in your local repo, you need to force the push. https://stackoverflow.com/questions/5509543/how-do-i-properly-force-a-git-push

[jcalbert]

OK I'm pretty sure I've got it now.

[jcalbert]

It's been a while but I think this still applies.

[mdhaber]

@steppi @bsdz @ragibson Ran across this and thought you might have an opinion. Shall I move this over to stats/_levy_stable/__init__.py? Where in the notes would it fit?

[mdhaber]

Miraculously, Sphinx rendered the table exactly as I would have hoped.

Here's a direct link to the artifact.
https://43418-1460385-gh.circle-artifacts.com/0/html-scipyorg/reference/generated/scipy.stats.levy_stable.html#scipy.stats.levy_stable

[bsdz]

Miraculously, Sphinx rendered the table exactly as I would have hoped.

Here's a direct link to the artifact. https://43418-1460385-gh.circle-artifacts.com/0/html-scipyorg/reference/generated/scipy.stats.levy_stable.html#scipy.stats.levy_stable

We do specialize for most of those cases. E.g.

scipy/scipy/stats/_levy_stable/__init__.py

Lines 190 to 211 in 5ec3c62

	if alpha == 2.0:
	# normal
	return _norm_pdf(x0 / np.sqrt(2)) / np.sqrt(2)
	elif alpha == 0.5 and beta == 1.0:
	# levy
	# since S(1/2, 1, gamma, delta; <x>) ==
	# S(1/2, 1, gamma, gamma + delta; <x0>).
	_x = x0 + 1
	if _x <= 0:
	return 0
	return 1 / np.sqrt(2 * np.pi * _x) / _x * np.exp(-1 / (2 * _x))
	elif alpha == 0.5 and beta == 0.0 and x0 != 0:
	# analytical solution [HO]
	S, C = sc.fresnel([1 / np.sqrt(2 * np.pi * np.abs(x0))])
	arg = 1 / (4 * np.abs(x0))
	return (
	np.sin(arg) * (0.5 - S[0]) + np.cos(arg) * (0.5 - C[0])
	) / np.sqrt(2 * np.pi * np.abs(x0) ** 3)
	elif alpha == 1.0 and beta == 0.0:
	# cauchy
	return 1 / (1 + x0 ** 2) / np.pi

That said, I feel it wise to recommend the equivalent distributions where possible. There might optimisations, better consideration of numerics etc. That table makes it clear.

[mdhaber]

Thanks. If that looks like a good place for it, I'll merge.

[ragibson]

I also agree that this information is useful and could be added to the documentation.

Wikipedia has a few more special cases (Landau and Holtsmark distributions) and analytic cases (using Bessel functions, hypergeometric functions, Fresnel integrals, etc.) for any particularly interested readers.
https://en.wikipedia.org/wiki/Stable_distribution#Special_cases
https://en.wikipedia.org/wiki/Stable_distribution#Other_analytic_cases

[ragibson]

I also think the current location in the notes is reasonable -- it essentially comes at the end of the distribution overview and right before the more technical implementation details.

[steppi]

This looks good to me. The information is useful and I also agree with its placement in the notes.

[mdhaber]

OK. I'll go ahead and merge as-is, and if we want to add additional information later, we can. Thanks @jcalbert and all!