Skip to content

Count trials in geometric, as the docstring and NumPy do - #196

Open
vlad-perevezentsev wants to merge 7 commits into
masterfrom
fix/geometric-numpy-convention
Open

Count trials in geometric, as the docstring and NumPy do#196
vlad-perevezentsev wants to merge 7 commits into
masterfrom
fix/geometric-numpy-convention

Conversation

@vlad-perevezentsev

@vlad-perevezentsev vlad-perevezentsev commented Oct 16, 2026 •
edited
Loading

Copy link
Copy Markdown
Collaborator

viRngGeometric returns the number of failures before the first success, with support 0, 1, 2, ... and mean 1 / p - 1, and that value was returned unchanged. The docstring documents the NumPy convention instead: "supported on the positive integers, k = 1, 2, ...", a PMF of (1 - p)^(k - 1) * p, and an example answering "How many trials succeeded after a single run?" with p.

geometric(1.0) returning 0 where NumPy returns 1 is just the most visible symptom: every p was off by one. Adding the successful trial fixes the whole distribution, and the special case for p == 1 now fills ones.

Closes #182

Changed

  • geometric samples are one larger than before for every p, so results for a fixed seed change accordingly. All other streams are unaffected.
  • The regression test records the new values; the added test pins the convention (support, P(1), mean) rather than a stream.

Notes

  • Values that would overflow int are left as they are — MKL already writes INT_MIN once 1 / p exceeds the range of its int output.

@vlad-perevezentsev
Count trials in geometric, as the docstring and NumPy do
eb69ce0 
viRngGeometric returns the number of failures, with support 0, 1, 2, ...
and mean 1 / p - 1, and that value was returned unchanged. The docstring
documents the NumPy convention instead: "supported on the positive
integers, k = 1, 2, ...", a PMF of (1 - p)^(k - 1) * p, and an example
answering "How many trials succeeded after a single run?" with p itself.
Measured over 200k samples at p = 0.35, the frequency of 1 was p * (1 - p)
= 0.2287, and p = 0.3487 was the frequency of 0.

This is not an edge case of p == 1, which merely made the mismatch obvious
by returning 0 where NumPy returns 1: every p was off by one. Adding the
successful trial fixes the whole distribution, so the special case for
p == 1 now fills ones. Values that would overflow int are left as they are;
they are already meaningless, as MKL writes INT_MIN once 1 / p exceeds the
range of its int output.

Samples for a given seed therefore grow by one, which the regression test
records. The new test pins the convention itself rather than a stream.
@vlad-perevezentsev vlad-perevezentsev self-assigned this Oct 16, 2026
@vlad-perevezentsev
Add gh-196 to changelog
b3b8f7a
@vlad-perevezentsev
Link the changelog entries to this PR
643d2a5 
gh-194 is the unrelated np.long issue. The entries were written before the
PR existed and kept the guessed number, the same slip as in gh-195.
antonwolfy
antonwolfy reviewed Oct 16, 2026
View reviewed changes

@antonwolfy antonwolfy left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Thanks, the fix looks right to me. I built it and compared against NumPy for p from 1.0 down to 1e-6, with scalar and array p, broadcasting and every brng: min, P(1), mean and variance all match now, and the full suite passes.

Could you add "Closes #182" to the description? That's the issue reporting this exact bug. Also, #195 touches the same changelog section and test_random.py, so whichever of the two lands second will need a rebase.

Comment thread CHANGELOG.md Outdated
Comment thread CHANGELOG.md Outdated
Comment thread mkl_random/src/mkl_distributions.cpp
Comment thread mkl_random/src/mkl_distributions.cpp Outdated
Comment thread mkl_random/tests/test_random.py Outdated
@antonwolfy antonwolfy added this to the 1.6.0 release milestone Oct 16, 2026
@vlad-perevezentsev
Keep MKL's overflow marker in geometric
2d7ce49 
Per review: the guard never fired. MKL writes INT_MIN on overflow, and
INT_MIN < INT_MAX, so the +1 turned it into INT_MIN + 1, which contradicts
the comment. Measured over 500k draws at p from 1e-9 down to 1e-12, the
only negative value present was INT_MIN + 1 and INT_MAX never appeared;
with the sign check the marker now survives unchanged.

The comment also explains why total and out exist: the chunk loop advances
res and decrements len, so without them the +1 pass would only cover the
last chunk for outputs above MKL_INT_MAX, which nothing tests.
@vlad-perevezentsev
Compare dtypes in the geometric ones check
4c45562 
Per review: assert_array_equal ignores dtypes by default, so the int32
expectation was not actually checked and an int64 or float result would
have passed. strict=True also pins the shape and needs NumPy 1.24, below
the 1.26.4 minimum.
@vlad-perevezentsev
Point the changelog entry at the old values
05fe07f 
Per review: "with a fixed seed they change accordingly" mostly repeated
"one larger", so the entry now says how to get the previous samples back.
vlad-perevezentsev added a commit that referenced this pull request Oct 16, 2026
@vlad-perevezentsev
Compare dtypes in the zeros and ones checks
473d1db 
Same point as the review on gh-196: assert_array_equal ignores dtypes by
default, so the int32 expectations were not actually checked and an int64
or float result would have passed. strict=True also pins the shape and
needs NumPy 1.24, below the 1.26.4 minimum.
antonwolfy
antonwolfy approved these changes Oct 16, 2026
View reviewed changes

@antonwolfy antonwolfy left a comment •
edited
Loading

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Thanks, this looks good now. I rebuilt 05fe07f and compared it against NumPy: chi-square against the exact PMF passes for p from 0.999 to 1e-5 across all brngs, scalar and array p and interfaces.numpy_random, p == 1 returns 1s in every form, and the output is exactly the old stream plus one, with no extra randomness consumed and MKL's INT_MIN overflow marker kept unchanged.

Two small things:

  • The description still doesn't have "Closes #182" (GitHub shows no linked issue), so #182 would stay open after the merge.
  • A correction to my overflow comment: the wrapped positive values I mentioned around p = 1e-10 don't exist. Re-measuring, every overflowed draw is exactly INT_MIN and the positive ones are genuine draws below 2^31, so the new comment is accurate as written.

Comment thread mkl_random/tests/test_random.py
@vlad-perevezentsev
Pin MKL's overflow marker in the geometric test
4858382 
Per review: the +1 pass skips negatives so the marker survives, but nothing
checked it. At p = 1e-300 every draw overflows, so the assertion is exact
rather than statistical; measured as all INT_MIN for p of 1e-100, 1e-300
and the smallest subnormal, on both the scalar and the array path.

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Reviewers

@antonwolfy antonwolfy antonwolfy approved these changes
@ndgrigorian ndgrigorian Awaiting requested review from ndgrigorian ndgrigorian is a code owner
@xaleryb xaleryb Awaiting requested review from xaleryb xaleryb is a code owner
@jharlow-intel jharlow-intel Awaiting requested review from jharlow-intel jharlow-intel is a code owner

Assignees

@vlad-perevezentsev vlad-perevezentsev

Labels

None yet

Projects

None yet

Milestone

1.6.0 release

Development

Successfully merging this pull request may close these issues.

geometric() returns failures before the first success (k >= 0) instead of trials (k >= 1)

2 participants

@vlad-perevezentsev @antonwolfy