'Returns an estimate of the scale parameter at the current
parameter value.'

def estimate_scale(self):

    if isinstance(self.family, (families.Binomial, families.Poisson, _Multinomial)):
        return 1.0
    endog = self.endog_li
    cached_means = self.cached_means
    nobs = self.nobs
    varfunc = self.family.variance
    scale = 0.0
    fsum = 0.0
    for i in range(self.num_group):
        if (len(endog[i]) == 0):
            continue
        (expval, _) = cached_means[i]
        f = (self.weights_li[i] if (self.weights is not None) else 1.0)
        sdev = np.sqrt(varfunc(expval))
        resid = ((endog[i] - expval) / sdev)
        scale += (f * np.sum((resid ** 2)))
        fsum += (f * len(endog[i]))
    scale /= ((fsum * (nobs - self.ddof_scale)) / float(nobs))
    return scale
