code: def _W(self, mu, weights, y=None):
          
          return sp.sparse.diags((self.link.gradient(mu, self.distribution)**2 *
                                  self.distribution.V(mu=mu) *
                                  weights ** -1)**-0.5)
docs: 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, optional
              does nothing. just for compatibility with ExpectileGAM
  
          Returns
          -------
          weights : sp..sparse array of shape (n_samples, n_samples)