docs: def resid_dev(self, endog, mu, scale=1.):
          r"""
          Binomial 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
  
          """
  
          mu = self.link._clean(mu)
          if np.shape(self.n) == () and self.n == 1:
              one = np.equal(endog, 1)
              return np.sign(endog-mu)*np.sqrt(-2 *
                                               np.log(one * mu + (1 - one) *
                                                      (1 - mu)))/scale
          else:
              return (np.sign(endog - mu) *
                      np.sqrt(2 * self.n *
                              (endog * np.log(endog/mu + 1e-200) +
                               (1 - endog) * np.log((1 - endog)/(1 - mu) + 1e-200)))/scale)