"""
The one parameter exponential family distributions used by GLM.
"""

# TODO: quasi, quasibinomial, quasipoisson
# see
# http://www.biostat.jhsph.edu/~qli/biostatistics_r_doc/library/stats/html/family.html
# for comparison to R, and McCullagh and Nelder


from statsmodels.compat.scipy import SP_LT_17

import inspect
import warnings

import numpy as np
from scipy import special, stats

from statsmodels.tools.sm_exceptions import (
    ValueWarning,
)

from . import links as L, varfuncs as V

FLOAT_EPS = np.finfo(float).eps


class Family:
    """
    The parent class for one-parameter exponential families.

    Parameters
    ----------
    link : a link function instance
        Link is the linear transformation function.
        See the individual families for available links.
    variance : a variance function
        Measures the variance as a function of the mean probabilities.
        See the individual families for the default variance function.
    check_link : bool
        If True (default), then and exception is raised if the link is invalid
        for the family.
        If False, then the link is not checked.

    See Also
    --------
    :ref:`links` : Further details on links.
    """

    # TODO: change these class attributes, use valid somewhere...
    valid = [-np.inf, np.inf]
    links = []

    def _setlink(self, link):
        """
        Helper method to set the link for a family.

        Raises a ``ValueError`` exception if the link is not available. Note
        that  the error message might not be that informative because it tells
        you that the link should be in the base class for the link function.

        See statsmodels.genmod.generalized_linear_model.GLM for a list of
        appropriate links for each family but note that not all of these are
        currently available.
        """
        # TODO: change the links class attribute in the families to hold
        # meaningful information instead of a list of links instances such as
        # [<statsmodels.family.links.Log object at 0x9a4240c>,
        #  <statsmodels.family.links.Power object at 0x9a423ec>,
        #  <statsmodels.family.links.Power object at 0x9a4236c>]
        # for Poisson...
        self._link = link
        if self._check_link:
            if not isinstance(link, L.Link):
                raise TypeError("The input should be a valid Link object.")
            if hasattr(self, "links"):
                validlink = max([isinstance(link, _) for _ in self.links])
                if not validlink:
                    msg = "Invalid link for family, should be in %s. (got %s)"
                    raise ValueError(msg % (repr(self.links), link))

    def _getlink(self):
        """
        Helper method to get the link for a family.
        """
        return self._link

    # link property for each family is a pointer to link instance
    link = property(_getlink, _setlink, doc="Link function for family")

    def __init__(self, link, variance, check_link=True):
        self._check_link = check_link
        if inspect.isclass(link):
            warnmssg = (
                "Calling Family(..) with a link class is not allowed. Use an "
                "instance of a link class instead."
            )
            raise TypeError(warnmssg)

        self.link = link
        self.variance = variance

    def starting_mu(self, y):
        r"""
        Starting value for mu in the IRLS algorithm.

        Parameters
        ----------
        y : ndarray
            The untransformed response variable.

        Returns
        -------
        mu_0 : ndarray
            The first guess on the transformed response variable.

        Notes
        -----
        .. math::

           \mu_0 = (Y + \overline{Y})/2

        Only the Binomial family takes a different initial value.
        """
        return (y + y.mean()) / 2.0

    def weights(self, mu):
        r"""
        Weights for IRLS steps

        Parameters
        ----------
        mu : array_like
            The transformed mean response variable in the exponential family

        Returns
        -------
        w : ndarray
            The weights for the IRLS steps

        Notes
        -----
        .. math::

           w = 1 / (g'(\mu)^2  * Var(\mu))
        """
        return 1.0 / (self.link.deriv(mu) ** 2 * self.variance(mu))

    def deviance(self, endog, mu, var_weights=1.0, freq_weights=1.0, scale=1.0):
        r"""
        The deviance function evaluated at (endog, mu, var_weights,
        freq_weights, scale) for the distribution.

        Deviance is usually defined as twice the loglikelihood ratio.

        Parameters
        ----------
        endog : array_like
            The endogenous response variable
        mu : array_like
            The inverse of the link function at the linear predicted values.
        var_weights : array_like
            1d array of variance (analytic) weights. The default is 1.
        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 : ndarray
            The value of deviance function defined below.

        Notes
        -----
        Deviance is defined

        .. math::

           D = 2\sum_i (freq\_weights_i * var\_weights *
           (llf(endog_i, endog_i) - llf(endog_i, \mu_i)))

        where y is the endogenous variable. The deviance functions are
        analytically defined for each family.

        Internally, we calculate deviance as:

        .. math::
           D = \sum_i freq\_weights_i * var\_weights * resid\_dev_i  / scale
        """
        resid_dev = self._resid_dev(endog, mu)
        return np.sum(resid_dev * freq_weights * var_weights / scale)

    def resid_dev(self, endog, mu, var_weights=1.0, scale=1.0):
        r"""
        The deviance residuals

        Parameters
        ----------
        endog : array_like
            The endogenous response variable
        mu : array_like
            The inverse of the link function at the linear predicted values.
        var_weights : array_like
            1d array of variance (analytic) weights. The default is 1.
        scale : float, optional
            An optional scale argument. The default is 1.

        Returns
        -------
        resid_dev : float
            Deviance residuals as defined below.

        Notes
        -----
        The deviance residuals are defined by the contribution D_i of
        observation i to the deviance as

        .. math::
           resid\_dev_i = sign(y_i-\mu_i) \sqrt{D_i}

        D_i is calculated from the _resid_dev method in each family.
        Distribution-specific documentation of the calculation is available
        there.
        """
        resid_dev = self._resid_dev(endog, mu)
        resid_dev *= var_weights / scale
        return np.sign(endog - mu) * np.sqrt(np.clip(resid_dev, 0.0, np.inf))

    def fitted(self, lin_pred):
        r"""
        Fitted values based on linear predictors lin_pred.

        Parameters
        ----------
        lin_pred : ndarray
            Values of the linear predictor of the model.
            :math:`X \cdot \beta` in a classical linear model.

        Returns
        -------
        mu : ndarray
            The mean response variables given by the inverse of the link
            function.
        """
        fits = self.link.inverse(lin_pred)
        return fits

    def predict(self, mu):
        """
        Linear predictors based on given mu values.

        Parameters
        ----------
        mu : ndarray
            The mean response variables

        Returns
        -------
        lin_pred : ndarray
            Linear predictors based on the mean response variables.  The value
            of the link function at the given mu.
        """
        return self.link(mu)

    def loglike_obs(self, endog, mu, var_weights=1.0, scale=1.0):
        r"""
        The log-likelihood function for each observation in terms of the fitted
        mean response for the distribution.

        Parameters
        ----------
        endog : ndarray
            Usually the endogenous response variable.
        mu : ndarray
            Usually but not always the fitted mean response variable.
        var_weights : array_like
            1d array of variance (analytic) weights. The default is 1.
        scale : float
            The scale parameter. The default is 1.

        Returns
        -------
        ll_i : float
            The value of the loglikelihood evaluated at
            (endog, mu, var_weights, scale) as defined below.

        Notes
        -----
        This is defined for each family. endog and mu are not restricted to
        ``endog`` and ``mu`` respectively.  For instance, you could call
        both ``loglike(endog, endog)`` and ``loglike(endog, mu)`` to get the
        log-likelihood ratio.
        """
        raise NotImplementedError

    def loglike(self, endog, mu, var_weights=1.0, freq_weights=1.0, scale=1.0):
        r"""
        The log-likelihood function in terms of the fitted mean response.

        Parameters
        ----------
        endog : ndarray
            Usually the endogenous response variable.
        mu : ndarray
            Usually but not always the fitted mean response variable.
        var_weights : array_like
            1d array of variance (analytic) weights. The default is 1.
        freq_weights : array_like
            1d array of frequency weights. The default is 1.
        scale : float
            The scale parameter. The default is 1.

        Returns
        -------
        ll : float
            The value of the loglikelihood evaluated at
            (endog, mu, var_weights, freq_weights, scale) as defined below.

        Notes
        -----
        Where :math:`ll_i` is the by-observation log-likelihood:

        .. math::
           ll = \sum(ll_i * freq\_weights_i)

        ``ll_i`` is defined for each family. endog and mu are not restricted
        to ``endog`` and ``mu`` respectively.  For instance, you could call
        both ``loglike(endog, endog)`` and ``loglike(endog, mu)`` to get the
        log-likelihood ratio.
        """
        ll_obs = self.loglike_obs(endog, mu, var_weights, scale)
        return np.sum(ll_obs * freq_weights)

    def resid_anscombe(self, endog, mu, var_weights=1.0, scale=1.0):
        r"""
        The Anscombe residuals

        Parameters
        ----------
        endog : ndarray
            The endogenous response variable
        mu : ndarray
            The inverse of the link function at the linear predicted values.
        var_weights : array_like
            1d array of variance (analytic) weights. The default is 1.
        scale : float, optional
            An optional argument to divide the residuals by sqrt(scale).
            The default is 1.

        See Also
        --------
        statsmodels.genmod.families.family.Family : `resid_anscombe` for the
          individual families for more information

        Notes
        -----
        Anscombe residuals are defined by

        .. math::
           resid\_anscombe_i = \frac{A(y)-A(\mu)}{A'(\mu)\sqrt{Var[\mu]}} *
           \sqrt(var\_weights)

        where :math:`A'(y)=v(y)^{-\frac{1}{3}}` and :math:`v(\mu)` is the
        variance function :math:`Var[y]=\frac{\phi}{w}v(mu)`.
        The transformation :math:`A(y)` makes the residuals more normal
        distributed.
        """
        raise NotImplementedError

    def _clean(self, x):
        """
        Helper function to trim the data so that it is in (0,inf)

        Notes
        -----
        The need for this function was discovered through usage and its
        possible that other families might need a check for validity of the
        domain.
        """
        return np.clip(x, FLOAT_EPS, np.inf)


class Poisson(Family):
    """
    Poisson exponential family.

    Parameters
    ----------
    link : a link instance, optional
        The default link for the Poisson family is the log link. Available
        links are log, identity, and sqrt. See statsmodels.families.links for
        more information.
    check_link : bool
        If True (default), then and exception is raised if the link is invalid
        for the family.
        If False, then the link is not checked.

    Attributes
    ----------
    Poisson.link : a link instance
        The link function of the Poisson instance.
    Poisson.variance : varfuncs instance
        ``variance`` is an instance of
        statsmodels.genmod.families.varfuncs.mu

    See Also
    --------
    statsmodels.genmod.families.family.Family : Parent class for all links.
    :ref:`links` : Further details on links.
    """

    links = [L.Log, L.Identity, L.Sqrt]
    variance = V.mu
    valid = [0, np.inf]
    safe_links = [
        L.Log,
    ]

    def __init__(self, link=None, check_link=True):
        if link is None:
            link = L.Log()
        super().__init__(link=link, variance=Poisson.variance, check_link=check_link)

    def _resid_dev(self, endog, mu):
        r"""
        Poisson deviance residuals

        Parameters
        ----------
        endog : ndarray
            The endogenous response variable.
        mu : ndarray
            The inverse of the link function at the linear predicted values.

        Returns
        -------
        resid_dev : float
            Deviance residuals as defined below.

        Notes
        -----
        .. math::

           resid\_dev_i = 2 * (endog_i * \ln(endog_i / \mu_i) -
           (endog_i - \mu_i))
        """
        endog_mu = self._clean(endog / mu)
        resid_dev = endog * np.log(endog_mu) - (endog - mu)
        return 2 * resid_dev

    def loglike_obs(self, endog, mu, var_weights=1.0, scale=1.0):
        r"""
        The log-likelihood function for each observation in terms of the fitted
        mean response for the Poisson distribution.

        Parameters
        ----------
        endog : ndarray
            Usually the endogenous response variable.
        mu : ndarray
            Usually but not always the fitted mean response variable.
        var_weights : array_like
            1d array of variance (analytic) weights. The default is 1.
        scale : float
            The scale parameter. The default is 1.

        Returns
        -------
        ll_i : float
            The value of the loglikelihood evaluated at
            (endog, mu, var_weights, scale) as defined below.

        Notes
        -----
        .. math::
            ll_i = var\_weights_i / scale * (endog_i * \ln(\mu_i) - \mu_i -
            \ln \Gamma(endog_i + 1))
        """
        return (
            var_weights / scale * (endog * np.log(mu) - mu - special.gammaln(endog + 1))
        )

    def resid_anscombe(self, endog, mu, var_weights=1.0, scale=1.0):
        r"""
        The Anscombe residuals

        Parameters
        ----------
        endog : ndarray
            The endogenous response variable
        mu : ndarray
            The inverse of the link function at the linear predicted values.
        var_weights : array_like
            1d array of variance (analytic) weights. The default is 1.
        scale : float, optional
            An optional argument to divide the residuals by sqrt(scale).
            The default is 1.

        Returns
        -------
        resid_anscombe : ndarray
            The Anscombe residuals for the Poisson family defined below

        Notes
        -----
        .. math::

           resid\_anscombe_i = (3/2) * (endog_i^{2/3} - \mu_i^{2/3}) /
           \mu_i^{1/6} * \sqrt(var\_weights)
        """
        resid = (
            (3 / 2.0)
            * (endog ** (2 / 3.0) - mu ** (2 / 3.0))
            / (mu ** (1 / 6.0) * scale**0.5)
        )
        resid *= np.sqrt(var_weights)
        return resid

    def get_distribution(self, mu, scale=1.0, var_weights=1.0):
        r"""
        Frozen Poisson distribution instance for given parameters

        Parameters
        ----------
        mu : ndarray
            Usually but not always the fitted mean response variable.
        scale : float
            The scale parameter is ignored.
        var_weights : array_like
            1d array of variance (analytic) weights. The default is 1.
            var_weights are ignored for Poisson.

        Returns
        -------
        distribution instance

        """

        return stats.poisson(mu)


class Gaussian(Family):
    """
    Gaussian exponential family distribution.

    Parameters
    ----------
    link : a link instance, optional
        The default link for the Gaussian family is the identity link.
        Available links are log, identity, and inverse.
        See statsmodels.genmod.families.links for more information.
    check_link : bool
        If True (default), then and exception is raised if the link is invalid
        for the family.
        If False, then the link is not checked.

    Attributes
    ----------
    Gaussian.link : a link instance
        The link function of the Gaussian instance
    Gaussian.variance : varfunc instance
        ``variance`` is an instance of
        statsmodels.genmod.families.varfuncs.constant

    See Also
    --------
    statsmodels.genmod.families.family.Family : Parent class for all links.
    :ref:`links` : Further details on links.
    """

    links = [L.Log, L.Identity, L.InversePower]
    variance = V.constant
    safe_links = links

    def __init__(self, link=None, check_link=True):
        if link is None:
            link = L.Identity()
        super().__init__(link=link, variance=Gaussian.variance, check_link=check_link)

    def _resid_dev(self, endog, mu):
        r"""
        Gaussian deviance residuals

        Parameters
        ----------
        endog : ndarray
            The endogenous response variable.
        mu : ndarray
            The inverse of the link function at the linear predicted values.

        Returns
        -------
        resid_dev : float
            Deviance residuals as defined below.

        Notes
        -----
        .. math::

           resid\_dev_i = (endog_i - \mu_i) ** 2
        """
        return (endog - mu) ** 2

    def loglike_obs(self, endog, mu, var_weights=1.0, scale=1.0):
        r"""
        The log-likelihood function for each observation in terms of the fitted
        mean response for the Gaussian distribution.

        Parameters
        ----------
        endog : ndarray
            Usually the endogenous response variable.
        mu : ndarray
            Usually but not always the fitted mean response variable.
        var_weights : array_like
            1d array of variance (analytic) weights. The default is 1.
        scale : float
            The scale parameter. The default is 1.

        Returns
        -------
        ll_i : float
            The value of the loglikelihood evaluated at
            (endog, mu, var_weights, scale) as defined below.

        Notes
        -----
        If the link is the identity link function then the
        loglikelihood function is the same as the classical OLS model.

        .. math::

           llf = -nobs / 2 * (\log(SSR) + (1 + \log(2 \pi / nobs)))

        where

        .. math::

           SSR = \sum_i (Y_i - g^{-1}(\mu_i))^2

        If the links is not the identity link then the loglikelihood
        function is defined as

        .. math::

           ll_i = -1 / 2 \sum_i  * var\_weights * ((Y_i - mu_i)^2 / scale +
                                                \log(2 * \pi * scale))
        """
        ll_obs = -var_weights * (endog - mu) ** 2 / scale
        ll_obs += -np.log(scale / var_weights) - np.log(2 * np.pi)
        ll_obs /= 2
        return ll_obs

    def resid_anscombe(self, endog, mu, var_weights=1.0, scale=1.0):
        r"""
        The Anscombe residuals

        Parameters
        ----------
        endog : ndarray
            The endogenous response variable
        mu : ndarray
            The inverse of the link function at the linear predicted values.
        var_weights : array_like
            1d array of variance (analytic) weights. The default is 1.
        scale : float, optional
            An optional argument to divide the residuals by sqrt(scale).
            The default is 1.

        Returns
        -------
        resid_anscombe : ndarray
            The Anscombe residuals for the Gaussian family defined below

        Notes
        -----
        For the Gaussian distribution, Anscombe residuals are the same as
        deviance residuals.

        .. math::

           resid\_anscombe_i = (Y_i - \mu_i) / \sqrt{scale} *
           \sqrt(var\_weights)
        """
        resid = (endog - mu) / scale**0.5
        resid *= np.sqrt(var_weights)
        return resid

    def get_distribution(self, mu, scale, var_weights=1.0):
        r"""
        Frozen Gaussian distribution instance for given parameters

        Parameters
        ----------
        mu : ndarray
            Usually but not always the fitted mean response variable.
        scale : float
            The scale parameter is required argument for get_distribution.
        var_weights : array_like
            1d array of variance (analytic) weights. The default is 1.

        Returns
        -------
        distribution instance

        """

        scale_n = scale / var_weights
        return stats.norm(loc=mu, scale=np.sqrt(scale_n))


class Gamma(Family):
    """
    Gamma exponential family distribution.

    Parameters
    ----------
    link : a link instance, optional
        The default link for the Gamma family is the inverse link.
        Available links are log, identity, and inverse.
        See statsmodels.genmod.families.links for more information.
    check_link : bool
        If True (default), then and exception is raised if the link is invalid
        for the family.
        If False, then the link is not checked.

    Attributes
    ----------
    Gamma.link : a link instance
        The link function of the Gamma instance
    Gamma.variance : varfunc instance
        ``variance`` is an instance of
        statsmodels.genmod.family.varfuncs.mu_squared

    See Also
    --------
    statsmodels.genmod.families.family.Family : Parent class for all links.
    :ref:`links` : Further details on links.
    """

    links = [L.Log, L.Identity, L.InversePower]
    variance = V.mu_squared
    safe_links = [
        L.Log,
    ]

    def __init__(self, link=None, check_link=True):
        if link is None:
            link = L.InversePower()
        super().__init__(link=link, variance=Gamma.variance, check_link=check_link)

    def _resid_dev(self, endog, mu):
        r"""
        Gamma deviance residuals

        Parameters
        ----------
        endog : ndarray
            The endogenous response variable.
        mu : ndarray
            The inverse of the link function at the linear predicted values.

        Returns
        -------
        resid_dev : float
            Deviance residuals as defined below.

        Notes
        -----
        .. math::

           resid\_dev_i = 2 * ((endog_i - \mu_i) / \mu_i -
           \log(endog_i / \mu_i))
        """
        endog_mu = self._clean(endog / mu)
        resid_dev = -np.log(endog_mu) + (endog - mu) / mu
        return 2 * resid_dev

    def loglike_obs(self, endog, mu, var_weights=1.0, scale=1.0):
        r"""
        The log-likelihood function for each observation in terms of the fitted
        mean response for the Gamma distribution.

        Parameters
        ----------
        endog : ndarray
            Usually the endogenous response variable.
        mu : ndarray
            Usually but not always the fitted mean response variable.
        var_weights : array_like
            1d array of variance (analytic) weights. The default is 1.
        scale : float
            The scale parameter. The default is 1.

        Returns
        -------
        ll_i : float
            The value of the loglikelihood evaluated at
            (endog, mu, var_weights, scale) as defined below.

        Notes
        -----
        .. math::

            ll_i = var\_weights_i / scale * (\ln(var\_weights_i * endog_i /
            (scale * \mu_i)) - (var\_weights_i * endog_i) /
            (scale * \mu_i)) - \ln \Gamma(var\_weights_i / scale) - \ln(\endog_i)

        Note on weights parameterization
        --------------------------------
        statsmodels follows the SPSS/SAS definition for variance weights:
        Var(endog_i) = scale * mu_i² / var_weights_i — i.e., the effective
        dispersion is scale / var_weights_i.

        McCullagh & Nelder (1989) treat var_weights_i as a pure multiplier
        on the log-likelihood with Var(endog_i) = scale * mu_i². This leads
        to a log-likelihood that differs only by a constant in
        (endog_i, var_weights_i, scale), not mu_i. So estimates of β and
        deviance are identical, but absolute log-likelihood values (e.g., AIC)
        will differ by a constant.

        References
        ----------
        McCullagh, P., and Nelder, J.A. (1989). *Generalized Linear Models*,
        2nd Edition. Chapman & Hall/CRC.
        """
        endog_mu = self._clean(endog / mu)
        weight_scale = var_weights / scale
        ll_obs = weight_scale * np.log(weight_scale * endog_mu)
        ll_obs -= weight_scale * endog_mu
        ll_obs -= special.gammaln(weight_scale) + np.log(endog)
        return ll_obs

        # in Stata scale is set to equal 1 for reporting llf
        # in R it's the dispersion, though there is a loss of precision vs.
        # our results due to an assumed difference in implementation

    def resid_anscombe(self, endog, mu, var_weights=1.0, scale=1.0):
        r"""
        The Anscombe residuals

        Parameters
        ----------
        endog : ndarray
            The endogenous response variable
        mu : ndarray
            The inverse of the link function at the linear predicted values.
        var_weights : array_like
            1d array of variance (analytic) weights. The default is 1.
        scale : float, optional
            An optional argument to divide the residuals by sqrt(scale).
            The default is 1.

        Returns
        -------
        resid_anscombe : ndarray
            The Anscombe residuals for the Gamma family defined below

        Notes
        -----
        .. math::

           resid\_anscombe_i = 3 * (endog_i^{1/3} - \mu_i^{1/3}) / \mu_i^{1/3}
           / \sqrt{scale} * \sqrt(var\_weights)
        """
        resid = (
            3 * (endog ** (1 / 3.0) - mu ** (1 / 3.0)) / mu ** (1 / 3.0) / scale**0.5
        )
        resid *= np.sqrt(var_weights)
        return resid

    def get_distribution(self, mu, scale, var_weights=1.0):
        r"""
        Frozen Gamma distribution instance for given parameters

        Parameters
        ----------
        mu : ndarray
            Usually but not always the fitted mean response variable.
        scale : float
            The scale parameter is required argument for get_distribution.
        var_weights : array_like
            1d array of variance (analytic) weights. The default is 1.

        Returns
        -------
        distribution instance

        """
        # combine var_weights with scale
        scale_ = scale / var_weights
        shape = 1 / scale_
        scale_g = mu * scale_
        return stats.gamma(shape, scale=scale_g)


class Binomial(Family):
    """
    Binomial exponential family distribution.

    Parameters
    ----------
    link : a link instance, optional
        The default link for the Binomial family is the logit link.
        Available links are logit, probit, cauchy, log, loglog, and cloglog.
        See statsmodels.genmod.families.links for more information.
    check_link : bool
        If True (default), then and exception is raised if the link is invalid
        for the family.
        If False, then the link is not checked.

    Attributes
    ----------
    Binomial.link : a link instance
        The link function of the Binomial instance
    Binomial.variance : varfunc instance
        ``variance`` is an instance of
        statsmodels.genmod.families.varfuncs.binary

    See Also
    --------
    statsmodels.genmod.families.family.Family : Parent class for all links.
    :ref:`links` : Further details on links.

    Notes
    -----
    endog for Binomial can be specified in one of three ways:
    A 1d array of 0 or 1 values, indicating failure or success
    respectively.
    A 2d array, with two columns. The first column represents the
    success count and the second column represents the failure
    count.
    A 1d array of proportions, indicating the proportion of
    successes, with parameter `var_weights` containing the
    number of trials for each row.
    """

    links = [
        L.Logit,
        L.Probit,
        L.Cauchy,
        L.Log,
        L.LogC,
        L.CLogLog,
        L.LogLog,
        L.Identity,
    ]
    variance = V.binary  # this is not used below in an effort to include n

    # Other safe links, e.g. cloglog and probit are subclasses
    safe_links = [L.Logit, L.CDFLink]

    def __init__(self, link=None, check_link=True):  # , n=1.):
        if link is None:
            link = L.Logit()
        # TODO: it *should* work for a constant n>1 actually, if freq_weights
        # is equal to n
        self.n = 1
        # overwritten by initialize if needed but always used to initialize
        # variance since endog is assumed/forced to be (0,1)
        super().__init__(
            link=link, variance=V.Binomial(n=self.n), check_link=check_link
        )

    def starting_mu(self, y):
        r"""
        The starting values for the IRLS algorithm for the Binomial family.
        A good choice for the binomial family is :math:`\mu_0 = (Y_i + 0.5)/2`
        """
        return (y + 0.5) / 2

    def initialize(self, endog, freq_weights):
        """
        Initialize the response variable.

        Parameters
        ----------
        endog : ndarray
            Endogenous response variable
        freq_weights : ndarray
            1d array of frequency weights

        Returns
        -------
        If `endog` is binary, returns `endog`

        If `endog` is a 2d array, then the input is assumed to be in the format
        (successes, failures) and
        successes/(success + failures) is returned.  And n is set to
        successes + failures.
        """
        # if not np.all(np.asarray(freq_weights) == 1):
        #     self.variance = V.Binomial(n=freq_weights)
        if endog.ndim > 1 and endog.shape[1] > 2:
            raise ValueError(
                "endog has more than 2 columns. The Binomial "
                "link supports either a single response variable "
                "or a paired response variable."
            )
        elif endog.ndim > 1 and endog.shape[1] > 1:
            y = endog[:, 0]
            # overwrite self.freq_weights for deviance below
            self.n = endog.sum(1)
            return y * 1.0 / self.n, self.n
        else:
            return endog, np.ones(endog.shape[0])

    def _resid_dev(self, endog, mu):
        r"""
        Binomial deviance residuals

        Parameters
        ----------
        endog : ndarray
            The endogenous response variable.
        mu : ndarray
            The inverse of the link function at the linear predicted values.

        Returns
        -------
        resid_dev : float
            Deviance residuals as defined below.

        Notes
        -----
        .. math::

           resid\_dev_i = 2 * n * (endog_i * \ln(endog_i /\mu_i) +
           (1 - endog_i) * \ln((1 - endog_i) / (1 - \mu_i)))
        """
        endog_mu = self._clean(endog / (mu + 1e-20))
        n_endog_mu = self._clean((1.0 - endog) / (1.0 - mu + 1e-20))
        resid_dev = endog * np.log(endog_mu) + (1 - endog) * np.log(n_endog_mu)
        return 2 * self.n * resid_dev

    def loglike_obs(self, endog, mu, var_weights=1.0, scale=1.0):
        r"""
        The log-likelihood function for each observation in terms of the fitted
        mean response for the Binomial distribution.

        Parameters
        ----------
        endog : ndarray
            Usually the endogenous response variable.
        mu : ndarray
            Usually but not always the fitted mean response variable.
        var_weights : array_like
            1d array of variance (analytic) weights. The default is 1.
