docs: def deviance(self, endog, mu, freq_weights=1, scale=1., axis=None):
          r'''
          Deviance function for either Bernoulli or Binomial data.
  
          Parameters
          ----------
          endog : array-like
              Endogenous response variable (already transformed to a probability
              if appropriate).
          mu : array
              Fitted mean response variable
          freq_weights : array-like
              1d array of frequency weights. The default is 1.
          scale : float, optional
              An optional scale argument. The default is 1.
  
          Returns
          --------
          deviance : float
              The deviance function as defined below
  
          '''
          if np.shape(self.n) == () and self.n == 1:
              one = np.equal(endog, 1)
              return -2 * np.sum((one * np.log(mu + 1e-200) + (1-one) *
                                 np.log(1 - mu + 1e-200)) * freq_weights, axis=axis)
  
          else:
              return 2 * np.sum(self.n * freq_weights *
                                (endog * np.log(endog/mu + 1e-200) +
                                 (1 - endog) * np.log((1 - endog) /
                                 (1 - mu) + 1e-200)), axis=axis)