Add Fisher noncentral hypergeometric distribution - #726
Conversation
| def _norm(good, bad, n, odds): | ||
| mode = FisherNoncentralHypergeometric.support_point(None, None, good, bad, n, odds) | ||
| mode_scale = FisherNoncentralHypergeometric._logweight(mode, good, bad, n, odds) | ||
| max_draws = pt.max(n) + 1 | ||
| arange = pt.arange(max_draws).reshape((max_draws, 1)) | ||
| log_weights = FisherNoncentralHypergeometric._logweight(arange, good, bad, n, odds) | ||
| scaled_log_sum = pt.logsumexp(log_weights - mode_scale, axis=0) | ||
| return mode_scale + scaled_log_sum | ||
|
|
||
| def logp(value, good, bad, n, odds): | ||
| return FisherNoncentralHypergeometric._logweight( | ||
| value, good, bad, n, odds | ||
| ) - FisherNoncentralHypergeometric._norm(good, bad, n, odds) | ||
|
|
||
| def _logweight(value, good, bad, n, odds): | ||
| fails_value = n - value | ||
| result = ( | ||
| logpow(odds, value) | ||
| - factln(value) | ||
| - factln(good - value) | ||
| - factln(fails_value) | ||
| - factln(bad - fails_value) | ||
| ) | ||
|
|
||
| lower_bound = n - bad | ||
| lower = pt.switch(pt.gt(lower_bound, 0), lower_bound, 0) | ||
| upper = pt.switch(pt.lt(good, n), good, n) | ||
| res = pt.switch( | ||
| pt.lt(value, lower), |
There was a problem hiding this comment.
why aren't these all in one function?
There was a problem hiding this comment.
_norm could be inlined and I am happy to do that if you prefer. _logweight needs to be called three separate times but I could make it a local inside logp.
| pt.switch( | ||
| pt.lt(value, n), | ||
| pt.logsumexp( | ||
| FisherNoncentralHypergeometric.logp(pt.arange(value + 1), good, bad, n, odds), |
There was a problem hiding this comment.
Isn't this potentially huge? Should it use a Scan? Is there nothing better? What does scipy do?
There was a problem hiding this comment.
This whole implementation roughly follows scipy which just sums over the PMF to get the CDF. (Note that _nchypergeom_gen only implements the PMF. There is no special implementation of the CDF which is just inherited from rv_discrete.)
Scipy's PMF calculation itself is taken from a C++ library that computes logs of the terms in the numerator subtracting a scale which is the log of the term corresponding to the mean to keep it numerically well behaved. This is the pattern I followed in this implementation with slight differences (I use the mode/support point.)
This line in particular is summing probabilities so it should not be huge. I am not sure if Scan offers some advantage here, but since pytensor has a builtin logsumexp, which is exactly what is needed here, I just reached for that.
ricardoV94
left a comment
ricardoV94
left a comment
There was a problem hiding this comment.
Looks good @kaylimekay
I left some minor notes.
2 participants
This PR adds the Fisher non-central hypergeometric distribution. It is modeled after the corresponding PR in the main pymc repo that adds the usual (central) hypergeometric distribution. It incorporates the scipy implementation as the random variable.
This addition was first proposed in a discussion on the main pymc repo.
This distribution is particularly useful in hit rate analyses when quantifying the odds that a given signal is able to correctly identify some notion of "winners" in data. The author has applied the code in this PR to problems of this type.