Iteratively re-weighted least squares estimation routine

    Parameters
    ----------
    y           : array
                  n*1, dependent variable

    x           : array
                  n*k, designs matrix of k independent variables

    family      : family object
                  probability models: Gaussian, Poisson, or Binomial

    offset      : array
                  n*1, the offset variable for each observation.

    y_fix       : array
                  n*1, the fixed intercept value of y for each observation

    ini_betas   : array
                  1*k, starting values for the k betas within the iteratively
                  weighted least squares routine

    tol         : float
                  tolerance for estimation convergence

    max_iter    : integer maximum number of iterations if convergence not met

    wi          : array
                  n*1, weights to transform observations from location i in GWR



    Returns
    -------

    betas       : array
                  k*1, estimated coefficients

    mu          : array
                  n*1, predicted y values

    wx          : array
                  n*1, final weights used for iwls for GLM

    n_iter      : integer
                  number of iterations that when iwls algorithm terminates

    w           : array
                  n*1, final weights used for iwls for GWR

    z           : array
                  iwls throughput

    v           : array
                  iwls throughput

    xtx_inv_xt  : array
                  iwls throughout to compute GWR hat matrix
                  [X'X]^-1 X'

def iwls(y, x, family, offset, y_fix,
         ini_betas=None, tol=1.0e-8, max_iter=200, wi=None):
    """
    Iteratively re-weighted least squares estimation routine

    Parameters
    ----------
    y           : array
                  n*1, dependent variable

    x           : array
                  n*k, designs matrix of k independent variables

    family      : family object
                  probability models: Gaussian, Poisson, or Binomial

    offset      : array
                  n*1, the offset variable for each observation.

    y_fix       : array
                  n*1, the fixed intercept value of y for each observation

    ini_betas   : array
                  1*k, starting values for the k betas within the iteratively
                  weighted least squares routine

    tol         : float
                  tolerance for estimation convergence

    max_iter    : integer maximum number of iterations if convergence not met

    wi          : array
                  n*1, weights to transform observations from location i in GWR



    Returns
    -------

    betas       : array
                  k*1, estimated coefficients

    mu          : array
                  n*1, predicted y values

    wx          : array
                  n*1, final weights used for iwls for GLM

    n_iter      : integer
                  number of iterations that when iwls algorithm terminates

    w           : array
                  n*1, final weights used for iwls for GWR

    z           : array
                  iwls throughput

    v           : array
                  iwls throughput

    xtx_inv_xt  : array
                  iwls throughout to compute GWR hat matrix
                  [X'X]^-1 X'

    """
    n_iter = 0
    diff = 1.0e6

    if ini_betas is None:
        betas = np.zeros((x.shape[1], 1), np.float)
    else:
        betas = ini_betas

    if isinstance(family, Binomial):
        y = family.link._clean(y)
    if isinstance(family, Poisson):
        y_off = y / offset
        y_off = family.starting_mu(y_off)
        v = family.predict(y_off)
        mu = family.starting_mu(y)
    else:
        mu = family.starting_mu(y)
        v = family.predict(mu)

    while diff > tol and n_iter < max_iter:
        n_iter += 1
        w = family.weights(mu)
        z = v + (family.link.deriv(mu) * (y - mu))
        w = np.sqrt(w)
        if not isinstance(x, np.ndarray):
            w = sp.csr_matrix(w)
            z = sp.csr_matrix(z)
        wx = spmultiply(x, w, array_out=False)
        wz = spmultiply(z, w, array_out=False)
        if wi is None:
            n_betas = _compute_betas(wz, wx)
        else:
            n_betas, xtx_inv_xt = _compute_betas_gwr(wz, wx, wi)
        v = spdot(x, n_betas)
        mu = family.fitted(v)

        if isinstance(family, Poisson):
            mu = mu * offset

        diff = min(abs(n_betas - betas))
        betas = n_betas

    if wi is None:
        return betas, mu, wx, n_iter
    else:
        return betas, mu, v, w, z, xtx_inv_xt, n_iter