'Tweedie Deviance Residual
Parameters
endog : array-like
`endog` is the response variable
mu : array-like
`mu` is the fitted value of the model
scale : float, optional
An optional argument to divide the residuals by sqrt(scale).
The default is 1.
Returns
resid_dev : array
The array of deviance residuals
Notes
When :math:`p = 1`,
.. math::
dev_i = \mu_i
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)
The deviance residual is then
.. math::
resid\_dev_i = sign(Y_i-\mu_i) * \sqrt{2 * dev_i / scale}'

def resid_dev(self, endog, mu, 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 = self._clean(endog)
        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.sign((endog - mu)) * np.sqrt(((2 * dev) / scale)))
