code: def _initial_estimate(self, y, modelmat):
          
  
          # do a simple initialization for LinearGAMs
          if isinstance(self, LinearGAM):
              n, m = modelmat.shape
              return np.ones(m) * np.sqrt(EPS)
  
          # transform the problem to the linear scale
          y = deepcopy(y).astype('float64')
          y[y == 0] += .01 # edge case for log link, inverse link, and logit link
          y[y == 1] -= .01 # edge case for logit link
  
          y_ = self.link.link(y, self.distribution)
          y_ = make_2d(y_, verbose=False)
          assert np.isfinite(y_).all(), "transformed response values should be well-behaved."
  
          # solve the linear problem
          return np.linalg.solve(load_diagonal(modelmat.T.dot(modelmat).A),
                                 modelmat.T.dot(y_))
docs: Makes an inital estimate for the model coefficients.
  
          For a LinearGAM we simply initialize to small coefficients.
  
          For other GAMs we transform the problem to the linear space
          and solve an unpenalized version.
  
          Parameters
          ---------
          y : array-like of shape (n,)
              containing target data
          modelmat : sparse matrix of shape (n, m)
              containing model matrix of the spline basis
  
          Returns
          -------
          coef : array of shape (m,) containing the initial estimate for the model
              coefficients
  
          Notes
          -----
              This method implements the suggestions in
              Wood, section 2.2.2 Geometry and IRLS convergence, pg 80