def resid_dev(self, endog, mu, scale=1.):
        """
        Gaussian deviance residuals

        Parameters
        -----------
        endog : array-like
            Endogenous response variable
        mu : array-like
            Fitted mean response variable
        scale : float, optional
            An optional argument to divide the residuals by scale. The default
            is 1.

        Returns
        -------
        resid_dev : array
            Deviance residuals as defined below

        """

        return (endog - mu) / np.sqrt(self.variance(mu)) / scale