Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

ENH: stats: Crystalball distribution #6798

Merged
merged 4 commits into from
Jul 30, 2017
Merged

Conversation

thomaskeck
Copy link
Contributor

The crystalball function is used in high-energy physics to parametrize data.
It is basically a Gaussian with a one-sided power-law tail.
https://en.wikipedia.org/wiki/Crystal_Ball_function

This is a split of, of another PR which contained several other distributions #6795

@ev-br
Copy link
Member

ev-br commented Nov 19, 2016

One issue with tests which fail is that higher moments do not exist for n=3.
It'd be best to calculate the expressions for the moments analytically, and add a _munp method to return either a correct value or an inf. Failing that, tests of higher moments need to be skipped for this distribution.

@thomaskeck
Copy link
Contributor Author

I'm working on the _munp implementation,
but it is more complex than I expected.

Calculating the formulas using mathematica was easy, but there are numerical instabilities and odd definitions of the functions which are required to calculate the n-th moment (e.g. hyp2f1).

So if munp is required I will need some more time (maybe I find some time next weekend to dive deeper into the mathematics behind hyp2f1).

@ev-br
Copy link
Member

ev-br commented Nov 27, 2016

From a 10000 feet: Mathematica likes to express things in terms of hyp2f1, even if there's a simpler expression available (or at least it did in version 3 or something). Moments can be left for later I think.

@ev-br ev-br added the needs-work Items that are pending response from the author label Jan 25, 2017
@thomaskeck
Copy link
Contributor Author

Fixed merge conflicts.
Implemented _munp.

@ev-br
You were correct the hyp2f1 function were unnecessary, I calculated the integral by hand
using some substitutions and the binomial expansion; the unittests check the calculated
values from _munp against the values calculated by mathematica (wolframalpha.com).


%(example)s
"""
def _pdf(self, x, alpha, n):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe better use some other variable name, #5982

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I renamed alpha to beta and n to m.

N = 1.0 / (n/alpha / (n-1) * np.exp(-alpha**2 / 2.0) + _norm_pdf_C * _norm_cdf(alpha))
return N * _lazywhere(np.atleast_1d(x > -alpha), (x, alpha, n),
f = lambda x, alpha, n: np.exp(-x**2 / 2),
f2 = lambda x, alpha, n: (n/alpha)**n * np.exp(-alpha**2 / 2.0) * (n/alpha - alpha - x)**(-n))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Line way too long, here and elsewhere.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I refactored the lazywhere and introduced a right-hand-side (rhs) and left-hand-side (lhs)
lambda function which is passed to f and f2 of _lazywhere.
The lines are now shorter.

@thomaskeck
Copy link
Contributor Author

Ping.

I fixed the latest merge conflicts, are there further suggestions or comments to this PR?

@thomaskeck
Copy link
Contributor Author

Pong?

@thomaskeck
Copy link
Contributor Author

Hello, this is my bi-monthly reminder for this PR.
I understand that the PR is probably forgotten and not shown anymore on the first page of open PRs.
But still I try to get some attention by this comment :-)
@ev-br Since you are provided me feedback in the past I mention you directly, which I should have done maybe earlier :-)

@pv pv merged commit 0b7932d into scipy:master Jul 30, 2017
@pv
Copy link
Member

pv commented Jul 30, 2017

Looks good to me now. Thanks, merging!

@pv pv removed the needs-work Items that are pending response from the author label Jul 30, 2017
@pv
Copy link
Member

pv commented Jul 30, 2017

Some issues remaining: gh-7687

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement A new feature or improvement scipy.stats
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants