Add HyperGeometric Distribution to pymc3.distributions.discrete #4108 #4249

Harivallabha · 2020-11-23T21:42:01Z

Creating a separate PR and linking it to #4108 because I had a hard time rebasing the old commits on top of master, too many conflicts. So rather than rebase, I just pulled latest master and again added these commits on top of it.

Included suggestions by @tirthasheshpatel. Thank you! :D
Removed _repr_latex_ implementations.

@Spaak

pre-commit failed

codecov · 2020-11-24T10:29:58Z

Codecov Report

Merging #4249 (3e9c589) into master (6ef0b2d) will increase coverage by 0.00%.
The diff coverage is 92.00%.

@@           Coverage Diff           @@
##           master    #4249   +/-   ##
=======================================
  Coverage   87.59%   87.60%           
=======================================
  Files          88       88           
  Lines       14316    14341   +25     
=======================================
+ Hits        12540    12563   +23     
- Misses       1776     1778    +2

Impacted Files	Coverage Δ
pymc3/distributions/discrete.py	`96.78% <91.66%> (-0.43%)`	⬇️
pymc3/distributions/__init__.py	`100.00% <100.00%> (ø)`

MarcoGorelli

Thanks @Harivallabha

Pre-commit failures can be fixed with

pre-commit run --files pymc3/distributions/__init__.py pymc3/distributions/discrete.py # list other modified files here

, or by enabling pre-commit locally:

pre-commit install

so the checks run automatically when you commit

MarcoGorelli · 2020-11-24T17:58:17Z

pymc3/tests/test_distributions_random.py

@@ -737,6 +742,9 @@ def ref_rand(size, alpha, mu):
    def test_geometric(self):
        pymc3_random_discrete(pm.Geometric, {"p": Unit}, size=500, fails=50, ref_rand=nr.geometric)

+    def test_hypergeometric(self):
+        pymc3_random_discrete(pm.HyperGeometric, {"N": Nat, "n": Nat, "k": Nat}, size=500, fails=50, ref_rand=nr.hypergeometric)


np.random.hypergeometric takes arguments named differently than N, n, and k - you may need to define a helper ref_rand function here, as is done above in test_negative_binomial

MarcoGorelli · 2020-11-24T17:59:46Z

pymc3/distributions/discrete.py

+            - betaln(n - value + 1, bad - n + value + 1)
+            - betaln(tot + 1, 1)
+        )
+        lower = tt.switch(tt.gt(n - N + k, 0), n - N + k, 0)


I'm a bit of a theano noob, but would

Suggested change

lower = tt.switch(tt.gt(n - N + k, 0), n - N + k, 0)

lower = tt.clip(n - N + k, 0, n - N + k)

work?

No, that diverges from the value returned by scipy.stats.hypergeom; However removing the bound entirely, and just returning the result seems to match

Spaak · 2020-12-02T13:42:54Z

@Harivallabha great to see this moving forward; if I understand the code correctly we're almost there, right? Any thoughts on the failing tests?

(So far this PR is still listed to go into the 3.10 release, which I hope we'll be able to make soon. If you think it's better to wait, then let me know of course.)

twiecki · 2020-12-02T14:04:00Z

I think we can leave this as option for 3.10 but not block on it.

MarcoGorelli · 2020-12-02T14:05:42Z

Any thoughts on the failing tests?

Just for reference: the CI has been flaky recently (seems fixed now though!), but the error here seems related to the changes:

________________ TestScalarParameterSamples.test_hypergeometric ________________
[gw0] linux -- Python 3.6.11 /usr/share/miniconda/envs/testenv/bin/python

self = <pymc3.tests.test_distributions_random.TestScalarParameterSamples object at 0x7f701ddb5cf8>

    def test_hypergeometric(self):
        def ref_rand(size, N, k, n):
            return st.hypergeom.rvs(M=N, n=k, N=n, size=size)
    
        pymc3_random_discrete(
            pm.HyperGeometric,
            {"N": Nat, "k": Nat, "n": Nat},
            size=500,
            fails=100,
>           ref_rand=ref_rand,
        )

pymc3/tests/test_distributions_random.py:754: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

dist = <class 'pymc3.distributions.discrete.HyperGeometric'>
paramdomains = {'N': <pymc3.tests.test_distributions.Domain object at 0x7f708b011e10>, 'k': <pymc3.tests.test_distributions.Domain object at 0x7f708b011e10>, 'n': <pymc3.tests.test_distributions.Domain object at 0x7f708b011e10>}
valuedomain = <pymc3.tests.test_distributions.Domain object at 0x7f708acfbe10>
ref_rand = <function TestScalarParameterSamples.test_hypergeometric.<locals>.ref_rand at 0x7f701ddf3e18>
size = 500, alpha = 0.05, fails = 100

    def pymc3_random_discrete(
        dist, paramdomains, valuedomain=Domain([0]), ref_rand=None, size=100000, alpha=0.05, fails=20
    ):
        model = build_model(dist, valuedomain, paramdomains)
        domains = paramdomains.copy()
        for pt in product(domains, n_samples=100):
            pt = pm.Point(pt, model=model)
            p = alpha
            # Allow Chisq test to fail (i.e., the samples be different)
            # a certain number of times.
            f = fails
            while p <= alpha and f > 0:
                o = model.named_vars["value"].random(size=size, point=pt)
                e = ref_rand(size=size, **pt)
                o = np.atleast_1d(o).flatten()
                e = np.atleast_1d(e).flatten()
                observed = dict(zip(*np.unique(o, return_counts=True)))
                expected = dict(zip(*np.unique(e, return_counts=True)))
                for e in expected.keys():
                    expected[e] = (observed.get(e, 0), expected[e])
                k = np.array([v for v in expected.values()])
                if np.all(k[:, 0] == k[:, 1]):
                    p = 1.0
                else:
                    _, p = st.chisquare(k[:, 0], k[:, 1])
                f -= 1
>           assert p > alpha, str(pt)
E           AssertionError: {'N': array(1), 'k': array(1), 'n': array(1)}
E           assert nan > 0.05

pymc3/tests/test_distributions_random.py:111: AssertionError

…discrete diverges from expected value

…st_distributions_random.py

…test passing

Harivallabha · 2020-12-03T13:00:17Z

Yes, please don't block on this.

So the match_with_scipy tests are passing but the random tests are still failing.

The reason for the failing random tests is that scipy.stats.hypergeom.rvs(M=N, n=k, N=n, size=size) returns a ValueError when N>M. So, for instance, when N, k, n = draw_values([self.N, self.k, self.n], point=point, size=size) returns [1, 1, 2], we have M < N (1 < 2), which results in a ValueError when the scipy rvs is called.

When I handle this, I keep getting the AssertionError (nan > 0.05) in pymc3_random_discrete that Marco has pointed out above. I'll push the latest code so that you can reproduce the error, but the random test is still failing, and I'm not sure why I'm getting the above AssertionError.

tirthasheshpatel · 2020-12-03T13:50:44Z

pymc3/tests/test_distributions_random.py

+        pymc3_random_discrete(
+            pm.HyperGeometric,
+            {"N": Nat, "k": Nat, "n": Nat},
+            size=500,
+            fails=50,
+            ref_rand=ref_rand,
+        )


Suggested change

pymc3_random_discrete(

pm.HyperGeometric,

{"N": Nat, "k": Nat, "n": Nat},

size=500,

fails=50,

ref_rand=ref_rand,

)

pymc3_random_discrete(

pm.HyperGeometric,

{"N": Domain([10, 11, 12, 13], "int64"),

"k": Domain([4, 5, 6, 7], "int64"),

"n": Domain([6, 7, 8, 9], "int64")},

size=500,

fails=50,

ref_rand=ref_rand,

)

You can create your own domain like this.

This should avoid the error that you are facing

Specifying domain fixed it. Wonderful, thanks!

tirthasheshpatel · 2020-12-03T13:57:10Z

pymc3/distributions/discrete.py

+
+    ========  =============================
+
+    Support   :math:`x in [max(0, n - \mathbb{N} + k), min(k, n)]`


Suggested change

Support :math:`x in [max(0, n - \mathbb{N} + k), min(k, n)]`

Support :math:`x \in \left[\max(0, n - N + k), \min(k, n)\right]`

Nitpick for proper latex rendering.

tirthasheshpatel

A very small final nitpick. Otherwise, it looks good to me!

pymc3/distributions/discrete.py

Remove stray newline Co-authored-by: Tirth Patel <[email protected]>

Spaak

Looks good to me!

Spaak · 2020-12-03T17:31:09Z

Thanks @Harivallabha and @tirthasheshpatel! 🎉 (Edit: and @MarcoGorelli of course! :) )

twiecki · 2020-12-03T18:20:07Z

Great job everyone!

twiecki previously approved these changes Nov 24, 2020

View reviewed changes

twiecki assigned MarcoGorelli Nov 24, 2020

Spaak added this to the 3.10 milestone Nov 24, 2020

MarcoGorelli suggested changes Nov 24, 2020

View reviewed changes

Harivallabha force-pushed the Add-HyperGeometric-Dist-2 branch from fc8351c to 9882ff6 Compare November 25, 2020 20:50

MarcoGorelli self-requested a review November 26, 2020 10:11

twiecki removed this from the 3.10 milestone Dec 2, 2020

Harivallabha added 7 commits December 3, 2020 18:19

Add HyperGeometric distribution to discrete.py; Add tests

ddf5617

Add HyperGeo to distirbutions/__init__.py

dfb0440

Fix minor linting issue

3426e36

Add ref_rand helper function. Clip lower in logp

ae9fb7c

Fix bug. Now pymc3_matches_scipy runs without error but pymc3_random_…

164a8f2

…discrete diverges from expected value

passes match with scipy test in test_distributions.py but fails in te…

f4edd8d

…st_distributions_random.py

Modify HyperGeom.random; Random test still failing. match_with_scipy …

6cd35fd

…test passing

Harivallabha force-pushed the Add-HyperGeometric-Dist-2 branch from 9882ff6 to 6cd35fd Compare December 3, 2020 13:05

rm stray print

71508a0

tirthasheshpatel reviewed Dec 3, 2020

View reviewed changes

Fix failing random test by specifying domain

b0da512

tirthasheshpatel approved these changes Dec 3, 2020

View reviewed changes

pymc3/distributions/discrete.py Outdated Show resolved Hide resolved

Harivallabha and others added 2 commits December 3, 2020 21:24

Update pymc3/distributions/discrete.py

a9f13bc

Remove stray newline Co-authored-by: Tirth Patel <[email protected]>

Add note in RELEASE-NOTES.md

3e9c589

Spaak approved these changes Dec 3, 2020

View reviewed changes

Spaak added this to the 3.10 milestone Dec 3, 2020

Spaak merged commit 198d13e into pymc-devs:master Dec 3, 2020

ricardoV94 mentioned this pull request Dec 21, 2020

HyperGeometric distribution gives wrong results #4366

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add HyperGeometric Distribution to pymc3.distributions.discrete #4108 #4249

Add HyperGeometric Distribution to pymc3.distributions.discrete #4108 #4249

Harivallabha commented Nov 23, 2020

codecov bot commented Nov 24, 2020 •

edited

Loading

MarcoGorelli left a comment

MarcoGorelli Nov 24, 2020

Harivallabha Nov 25, 2020

MarcoGorelli Nov 24, 2020

Harivallabha Nov 25, 2020

Spaak commented Dec 2, 2020

twiecki commented Dec 2, 2020

MarcoGorelli commented Dec 2, 2020

Harivallabha commented Dec 3, 2020

tirthasheshpatel Dec 3, 2020

tirthasheshpatel Dec 3, 2020

Harivallabha Dec 3, 2020

tirthasheshpatel Dec 3, 2020

Harivallabha Dec 3, 2020

tirthasheshpatel left a comment

Spaak left a comment

Spaak commented Dec 3, 2020 •

edited

Loading

twiecki commented Dec 3, 2020

	lower = tt.switch(tt.gt(n - N + k, 0), n - N + k, 0)
	lower = tt.clip(n - N + k, 0, n - N + k)


		======== =============================

		Support :math:`x in [max(0, n - \mathbb{N} + k), min(k, n)]`

	Support :math:`x in [max(0, n - \mathbb{N} + k), min(k, n)]`
	Support :math:`x \in \left[\max(0, n - N + k), \min(k, n)\right]`

Add HyperGeometric Distribution to pymc3.distributions.discrete #4108 #4249

Add HyperGeometric Distribution to pymc3.distributions.discrete #4108 #4249

Conversation

Harivallabha commented Nov 23, 2020

codecov bot commented Nov 24, 2020 • edited Loading

Codecov Report

MarcoGorelli left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Spaak commented Dec 2, 2020

twiecki commented Dec 2, 2020

MarcoGorelli commented Dec 2, 2020

Harivallabha commented Dec 3, 2020

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tirthasheshpatel left a comment

Choose a reason for hiding this comment

Spaak left a comment

Choose a reason for hiding this comment

Spaak commented Dec 3, 2020 • edited Loading

twiecki commented Dec 3, 2020

codecov bot commented Nov 24, 2020 •

edited

Loading

Spaak commented Dec 3, 2020 •

edited

Loading