'Calculate expected value of a function with respect to the distribution
for discrete distribution.
Parameters
func : callable, optional
Function for which the expectation value is calculated.
Takes only one argument.
The default is the identity mapping f(k) = k.
args : tuple, optional
Shape parameters of the distribution.
loc : float, optional
Location parameter.
Default is 0.
lb, ub : int, optional
Lower and upper bound for the summation, default is set to the
support of the distribution, inclusive (``ul <= k <= ub``).
conditional : bool, optional
If true then the expectation is corrected by the conditional
probability of the summation interval. The return value is the
expectation of the function, `func`, conditional on being in
the given interval (k such that ``ul <= k <= ub``).
Default is False.
maxcount : int, optional
Maximal number of terms to evaluate (to avoid an endless loop for
an infinite sum). Default is 1000.
tolerance : float, optional
Absolute tolerance for the summation. Default is 1e-10.
chunksize : int, optional
Iterate over the support of a distributions in chunks of this size.
Default is 32.
Returns
expect : float
Expected value.
Notes
For heavy-tailed distributions, the expected value may or may not exist,
depending on the function, `func`. If it does exist, but the sum converges
slowly, the accuracy of the result may be rather low. For instance, for
``zipf(4)``, accuracy for mean, variance in example is only 1e-5.
increasing `maxcount` and/or `chunksize` may improve the result, but may also
make zipf very slow.
The function is not vectorized.'

def expect(self, func=None, args=(), loc=0, lb=None, ub=None, conditional=False, maxcount=1000, tolerance=1e-10, chunksize=32):

    if (func is None):
        def fun(x):
            return ((x + loc) * self._pmf(x, *args))
    else:
        def fun(x):
            return (func((x + loc)) * self._pmf(x, *args))
    self._argcheck(*args)
    if (lb is None):
        lb = self.a
    else:
        lb = (lb - loc)
    if (ub is None):
        ub = self.b
    else:
        ub = (ub - loc)
    if conditional:
        invfac = (self.sf((lb - 1), *args) - self.sf(ub, *args))
    else:
        invfac = 1.0
    x0 = self.ppf(0.5, *args)
    res = _expect(fun, lb, ub, x0, self.inc, maxcount, tolerance, chunksize)
    return (res / invfac)
