Added Breit-Wigner distribution #6419
Conversation
|
A few issues I need to fix:
Expected: Got:
|
|
|
||
| _argcheck : bool | ||
| Whether shape parameters are legitimate | ||
| """ |
There was a problem hiding this comment.
Nitpick: There's no need to add a full-fledged docstring to a private, computational method. A simple comment would suffice IMO (and even then, the one-line code is self-documenting enough, I'd think)
|
First and foremost, my question about the use case (#6414) still stands. On a technical level, this is a very good start! As for doctests, use |
ev-br
added
scipy.stats
enhancement
| >>> [BW.moment(n) for n in range(1, 6)] | ||
| [122.11538219739913, 15674.920254744189, inf, 244138132.80994478, 80583978992.641495] | ||
|
|
||
| """ |
There was a problem hiding this comment.
Have a look at the docstrings of other distributions, they have magic substitution keys (compare to online docs, where the keys are expanded). I probably makes sense to reuse the standard format.
There was a problem hiding this comment.
The docstring also needs a See Also section with a reference to halfcauchy.
|
Any news on this? Would love this feature! |
|
The unanswered question is about whether And indeed, as a function of the square of the 4-momentum, E^2, Breit-Wigner formula is almost a Cauchy distribution, the only difference being that the support of the distribution is [0, \inf], so it's a halfcauchy. Can you sketch a use case where a distribution of E^2 is not enough and you really need a distribution in terms of E itself? Once that's decided, this looks close to good to go. |
| >>> samples = BW.rvs(1000) | ||
| >>> hist_plot = plt.hist(samples, normed=True, bins=100, range=[0, 200]) | ||
|
|
||
| >>> plt.show() |
There was a problem hiding this comment.
Explicit example largely duplicates the autogenerated. Suggest deduplicating.
| The CDF was found by Mathematica: | ||
|
|
||
| pdf = k/((m^2 - mass^2)^2 + mass^4*alpha^2) | ||
| cdf = Integrate[pdf, m] |
There was a problem hiding this comment.
This comment belongs to the comment of _cdf, should not be in the main docstring.
| We require mass > 0 and width > 0. | ||
|
|
||
| The width -> 0 limit is well-defined mathematically - it's proportional | ||
| to a Dirac function. In the m -> 0 limit, the pdf vanishes. |
|
|
||
| The CDF was found by Mathematica: | ||
|
|
||
| pdf = k/((m^2 - mass^2)^2 + mass^4*alpha^2) |
There was a problem hiding this comment.
Most other distribution have TeX formulas here.
| arg_1 = complex(-1)**(1. / 4.) / (-1j + alpha_)**0.5 * m / mass | ||
| arg_2 = complex(-1)**(3. / 4.) / (1j + alpha_)**0.5 * m / mass | ||
|
|
||
| shape_ = -1j * np.arctan(arg_1) / (-1j + alpha_)**0.5 - np.arctan(arg_2) / (1j + alpha_)**0.5 |
There was a problem hiding this comment.
I'm pretty sure this formula could be simplified further to avoid complex arithmetics.
I think it is, see #6414 (comment). |
|
I don't have time for this, sorry, but would happy to see it make it into scipy :) |
I've added the relativistic Breit-Wigner distribution to
scipy.stats. The_pdfand_cdfare analytic; other attributes are numerical and provided by the inheritedrv_continuousclass.As for the other distributions, I define a
breit_wigner_genclass, and then abreit_wignerinstance of that class, which should be imported by a user.