'Returns the value of the deviance function.
Parameters
endog : array-like
Endogenous response variable
mu : array-like
Fitted mean response variable
iweights : array-like
1d array of weights. The default is 1.
scale : float, optional
An optional scale argument. The default is 1.
Returns
deviance : float
Deviance function as defined below
Notes
When :math:`p = 1`,
.. math::
dev_i = \mu
when :math:`Y_i = 0` and
.. math::
dev_i = Y_i * \log(Y_i / \mu_i) + (\mu_i - Y_i)
otherwise.
When :math:`p = 2`,
.. math::
dev_i =  (Y_i - \mu_i) / \mu_i - \log(Y_i / \mu_i)
For all other p,
.. math::
dev_i = Y_i^{2 - p} / ((1 - p) * (2 - p)) -
Y_i * \mu_i^{1 - p} / (1 - p) + \mu_i^{2 - p} /
(2 - p)
Once :math:`dev_i` is calculated, then deviance is calculated as
.. math::
D = \sum{2 * iweights * dev_i / scale}'

def deviance(self, endog, mu, iweights=1.0, scale=1.0):

    p = self.var_power
    if (p == 1):
        dev = np.where((endog == 0), mu, ((endog * np.log((endog / mu))) + (mu - endog)))
    elif (p == 2):
        endog1 = np.clip(endog, FLOAT_EPS, np.inf)
        dev = (((endog - mu) / mu) - np.log((endog1 / mu)))
    else:
        dev = ((((endog ** (2 - p)) / ((1 - p) * (2 - p))) - ((endog * (mu ** (1 - p))) / (1 - p))) + ((mu ** (2 - p)) / (2 - p)))
    return np.sum((((2 * iweights) * dev) / scale))
