code: def test_nchypergeom_wallenius_naive(self):
      # test against a very simple implementation
  
      np.random.seed(2)
      shape = (2, 4, 3)
      max_m = 100
      m1 = np.random.randint(1, max_m, size=shape)
      m2 = np.random.randint(1, max_m, size=shape)
      N = m1 + m2
      n = randint.rvs(0, N, size=N.shape)
      xl = np.maximum(0, n - m2)
      xu = np.minimum(n, m1)
      x = randint.rvs(xl, xu, size=xl.shape)
      w = np.random.rand(*x.shape) * 2
  
      def support(N, m1, n, w):
          m2 = N - m1
          xl = np.maximum(0, n - m2)
          xu = np.minimum(n, m1)
          return xl, xu
  
      @np.vectorize
      def mean(N, m1, n, w):
          m2 = N - m1
          xl, xu = support(N, m1, n, w)
  
          def fun(u):
              return u / m1 + (1 - (n - u) / m2) ** w - 1
  
          return root_scalar(fun, bracket=(xl, xu)).root
  
      assert_allclose(
          nchypergeom_wallenius.mean(N, m1, n, w),
          mean(N, m1, n, w),
          rtol=2e-2,
      )
  
      @np.vectorize
      def variance(N, m1, n, w):
          m2 = N - m1
          u = mean(N, m1, n, w)
          a = u * (m1 - u)
          b = (n - u) * (u + m2 - n)
          return N * a * b / ((N - 1) * (m1 * b + m2 * a))
  
      assert_allclose(
          nchypergeom_wallenius.stats(N, m1, n, w, moments="v"),
          variance(N, m1, n, w),
          rtol=5e-2,
      )
  
      @np.vectorize
      def pmf(x, N, m1, n, w):
          m2 = N - m1
          xl, xu = support(N, m1, n, w)
  
          def integrand(t):
              D = w * (m1 - x) + (m2 - (n - x))
              res = (1 - t ** (w / D)) ** x * (1 - t ** (1 / D)) ** (
                  n - x
              )
              return res
  
          def f(x):
              t1 = special_binom(m1, x)
              t2 = special_binom(m2, n - x)
              the_integral = quad(
                  integrand, 0, 1, epsrel=1e-16, epsabs=1e-16
              )
              return t1 * t2 * the_integral[0]
  
          return f(x)
  
      pmf0 = pmf(x, N, m1, n, w)
      pmf1 = nchypergeom_wallenius.pmf(x, N, m1, n, w)
  
      atol, rtol = 1e-6, 1e-6
      i = np.abs(pmf1 - pmf0) < atol + rtol * np.abs(pmf0)
      assert (
          i.sum() > np.prod(shape) / 2
      )  # works at least half the time
  
      # for those that fail, discredit the naive implementation
      for N, m1, n, w in zip(N[~i], m1[~i], n[~i], w[~i]):
          # get the support
          m2 = N - m1
          xl, xu = support(N, m1, n, w)
          x = np.arange(xl, xu + 1)
  
          # calculate sum of pmf over the support
          # the naive implementation is very wrong in these cases
          assert pmf(x, N, m1, n, w).sum() < 0.5
          assert_allclose(
              nchypergeom_wallenius.pmf(x, N, m1, n, w).sum(), 1
          )
  
imports: import pytest
  from scipy.stats import betabinom, hypergeom, nhypergeom, bernoulli, boltzmann, skellam, zipf, zipfian, binom, nbinom, nchypergeom_fisher, nchypergeom_wallenius, randint
  import numpy
  from numpy.testing import assert_almost_equal, assert_equal, assert_allclose
  from scipy.special import binom
  import pytest
  from scipy.optimize import root_scalar
  from scipy.integrate import quad
label_desc: Use image node_id 3 for calling the TestNCH obj's underlying member method code with example usage: obj.test_nchypergeom_wallenius_naive() and returns: xl,  xu, root_scalar, unknown, f, unknown, res, unknown