compute the PIRLS weights for model predictions.

        TODO lets verify the formula for this.
        if we use the square root of the mu with the stable opt,
        we get the same results as when we use non-sqrt mu with naive opt.

        this makes me think that they are equivalent.

        also, using non-sqrt mu with stable opt gives very small edofs for even lam=0.001
        and the parameter variance is huge. this seems strange to me.

        computed [V * d(link)/d(mu)] ^(-1/2) by hand and the math checks out as hoped.

        ive since moved the square to the naive pirls method to make the code modular.

        Parameters
        ---------
        mu : array-like of shape (n_samples,)
            expected value of the targets given the model and inputs
        weights : array-like of shape (n_samples,)
            containing sample weights
        y = array-like of shape (n_samples,) or None, default None
            useful for computing the asymmetric weight.

        Returns
        -------
        weights : scipy.sparse array of shape (n_samples, n_samples)

def _W(self, mu, weights, y=None):
        """
        compute the PIRLS weights for model predictions.

        TODO lets verify the formula for this.
        if we use the square root of the mu with the stable opt,
        we get the same results as when we use non-sqrt mu with naive opt.

        this makes me think that they are equivalent.

        also, using non-sqrt mu with stable opt gives very small edofs for even lam=0.001
        and the parameter variance is huge. this seems strange to me.

        computed [V * d(link)/d(mu)] ^(-1/2) by hand and the math checks out as hoped.

        ive since moved the square to the naive pirls method to make the code modular.

        Parameters
        ---------
        mu : array-like of shape (n_samples,)
            expected value of the targets given the model and inputs
        weights : array-like of shape (n_samples,)
            containing sample weights
        y = array-like of shape (n_samples,) or None, default None
            useful for computing the asymmetric weight.

        Returns
        -------
        weights : scipy.sparse array of shape (n_samples, n_samples)
        """
        # asymmetric weight
        asym = (y > mu) * self.expectile + (y <= mu) * (1 - self.expectile)

        return sp.sparse.diags((self.link.gradient(mu, self.distribution)**2 *
                                self.distribution.V(mu=mu) *
                                weights ** -1)**-0.5 * asym**0.5)