code: def test_betanbinom_kurtosis(self, n, a, b, ref):
      # reference values were computed via mpmath
      # from mpmath import mp
      # def kurtosis_betanegbinom(n, a, b):
      #     n = mp.mpf(n)
      #     a = mp.mpf(a)
      #     b = mp.mpf(b)
      #     four = mp.mpf(4.)
      #     mean = n * b / (a - mp.one)
      #     var = (n * b * (n + a - 1.) * (a + b - 1.)
      #            / ((a - 2.) * (a - 1.)**2.))
      #     def f(k):
      #         return (mp.binomial(n + k - mp.one, k)
      #                 * mp.beta(a + n, b + k) / mp.beta(a, b)
      #                 * (k - mean)**four)
      #      fourth_moment = mp.nsum(f, [0, mp.inf])
      #      return float(fourth_moment/var**2 - 3.)
      assert_allclose(
          betanbinom.stats(n, a, b, moments="k"), ref, rtol=3e-15
      )
  
imports: import pytest
  import itertools
  from scipy.stats import betabinom, betanbinom, 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, suppress_warnings
  from scipy.special import binom
  from scipy.optimize import root_scalar
  from scipy.integrate import quad
label_desc: Use image node_id 2 for calling the TestBetaNBinom obj's underlying member method code with example usage: obj.test_betanbinom_kurtosis(n, a, b, ref) without return types