code: def phi(self, y, mu, edof, weights):
          
          if self._known_scale:
              return self.scale
          else:
              return (np.sum(weights * self.V(mu)**-1 * (y - mu)**2) /
                      (len(mu) - edof))
docs: GLM scale parameter.
          for Binomial and Poisson families this is unity
          for Normal family this is variance
  
          Parameters
          ----------
          y : array-like of length n
              target values
          mu : array-like of length n
              expected values
          edof : float
              estimated degrees of freedom
          weights : array-like shape (n,) or None, default: None
              sample weights
              if None, defaults to array of ones
  
          Returns
          -------
          scale : estimated model scale