Skip to content

BUG: random: beta can generate nan when the parameters are extremely small. #24266

Closed

Description

@WarrenWeckesser

Describe the issue:

This issue first came up in #24210. I'm creating a new issue to focus on beta.

The beta method of random.Generator can generate nan when both parameters are extremely small. "Extremely small" means subnormal, or normal but only a small multiple the smallest normal float. For example, in the following sample of size 500000000, 413 nan values were generated.

In [1]: rng = np.random.default_rng()

In [2]: tiny = np.finfo(np.float64).tiny

In [3]: x = rng.beta(2*tiny, 1.5*tiny, size=500000000)

In [4]: np.isnan(x).sum()
Out[4]: np.int64(413)

The problem is in the C function random_beta. In this loop

while (1) {
U = next_double(bitgen_state);
V = next_double(bitgen_state);
X = pow(U, 1.0 / a);
Y = pow(V, 1.0 / b);
XpY = X + Y;
/* Reject if both U and V are 0.0, which is approx 1 in 10^106 */
if ((XpY <= 1.0) && (U + V > 0.0)) {
if (XpY > 0) {
return X / XpY;
} else {
double logX = log(U) / a;
double logY = log(V) / b;
double logM = logX > logY ? logX : logY;
logX -= logM;
logY -= logM;
return exp(logX - log(exp(logX) + exp(logY)));
}
}
}

when a and b are sufficiently small, the both variables logX and logY can be -inf. Then the subtraction -inf - (-inf) results in nan.

Reproduce the code example:

See above.

Error message:

No response

Runtime information:

In [5]: import sys, numpy; print(numpy.__version__); print(sys.version)
1.25.0rc1+630.gd50fc570a9
3.11.4 (main, Jul  3 2023, 14:49:40) [GCC 11.3.0]

In [6]: print(numpy.show_runtime())
[{'numpy_version': '1.25.0rc1+630.gd50fc570a9',
  'python': '3.11.4 (main, Jul  3 2023, 14:49:40) [GCC 11.3.0]',
  'uname': uname_result(system='Linux', node='pop-os', release='6.2.6-76060206-generic', version='#202303130630~1689015125~22.04~ab2190e SMP PREEMPT_DYNAMIC Mon J', machine='x86_64')},
 {'simd_extensions': {'baseline': [], 'found': [], 'not_found': []}},
 {'architecture': 'Zen',
  'filepath': '/usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so',
  'internal_api': 'openblas',
  'num_threads': 24,
  'prefix': 'libopenblas',
  'threading_layer': 'pthreads',
  'user_api': 'blas',
  'version': '0.3.20'}]
None

Context for the issue:

No response

Activity

  1. added 2 commits that reference this issue on Jul 26, 2023
    0f7c3c8
    a0b45c5
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions