'Return Pearson\'s X^2 estimate of scale.'

def estimate_scale(self, Y=None):

    if (Y is None):
        Y = self.Y
    resid = (Y - self.results.mu)
    return ((np.power(resid, 2) / self.family.variance(self.results.mu)).sum() / self.df_resid)
