Simulate from the posterior of the coefficients.

        NOTE: A `gridsearch` is done `n_bootstraps` many times, so keep
        `n_bootstraps` small. Make `n_bootstraps < n_draws` to take advantage
        of the expensive bootstrap samples of the smoothing parameters.

        Parameters
        -----------
        X : array of shape (n_samples, m_features)
              input data

        y : array of shape (n_samples,)
              response vector

        weights : np.array of shape (n_samples,)
            sample weights

        n_draws : positive int, optional (default=100
            The number of samples to draw from the posterior distribution of
            the coefficients and smoothing parameters

        n_bootstraps : positive int, optional (default=1
            The number of bootstrap samples to draw from simulations of the
            response (from the already fitted model) to estimate the
            distribution of the smoothing parameters given the response data.
            If `n_bootstraps` is 1, then only the already fitted model's
            smoothing parameters is used.

        objective : string, optional (default='auto'
            metric to optimize in grid search. must be in
            ['AIC', 'AICc', 'GCV', 'UBRE', 'auto']
            if 'auto', then grid search will optimize GCV for models with
            unknown scale and UBRE for models with known scale.

        Returns
        -------
        coef_samples : array of shape (n_draws, n_samples)
            Approximate simulations of the coefficients drawn from the
            posterior distribution of the coefficients and smoothing
            parameters given the response data

        References
        ----------
        Simon N. Wood, 2006. Generalized Additive Models: an introduction with
        R. Section 4.9.3 (pages 198–199) and Section 5.4.2 (page 256–257).

def _sample_coef(self, X, y, weights=None, n_draws=100, n_bootstraps=1,
                     objective='auto'):
        """Simulate from the posterior of the coefficients.

        NOTE: A `gridsearch` is done `n_bootstraps` many times, so keep
        `n_bootstraps` small. Make `n_bootstraps < n_draws` to take advantage
        of the expensive bootstrap samples of the smoothing parameters.

        Parameters
        -----------
        X : array of shape (n_samples, m_features)
              input data

        y : array of shape (n_samples,)
              response vector

        weights : np.array of shape (n_samples,)
            sample weights

        n_draws : positive int, optional (default=100
            The number of samples to draw from the posterior distribution of
            the coefficients and smoothing parameters

        n_bootstraps : positive int, optional (default=1
            The number of bootstrap samples to draw from simulations of the
            response (from the already fitted model) to estimate the
            distribution of the smoothing parameters given the response data.
            If `n_bootstraps` is 1, then only the already fitted model's
            smoothing parameters is used.

        objective : string, optional (default='auto'
            metric to optimize in grid search. must be in
            ['AIC', 'AICc', 'GCV', 'UBRE', 'auto']
            if 'auto', then grid search will optimize GCV for models with
            unknown scale and UBRE for models with known scale.

        Returns
        -------
        coef_samples : array of shape (n_draws, n_samples)
            Approximate simulations of the coefficients drawn from the
            posterior distribution of the coefficients and smoothing
            parameters given the response data

        References
        ----------
        Simon N. Wood, 2006. Generalized Additive Models: an introduction with
        R. Section 4.9.3 (pages 198–199) and Section 5.4.2 (page 256–257).
        """
        if not self._is_fitted:
            raise AttributeError('GAM has not been fitted. Call fit first.')
        if n_bootstraps < 1:
            raise ValueError('n_bootstraps must be >= 1;'
                             ' got {}'.format(n_bootstraps))
        if n_draws < 1:
            raise ValueError('n_draws must be >= 1;'
                             ' got {}'.format(n_draws))

        coef_bootstraps, cov_bootstraps = (
            self._bootstrap_samples_of_smoothing(X, y, weights=weights,
                                                 n_bootstraps=n_bootstraps,
                                                 objective=objective))
        coef_draws = self._simulate_coef_from_bootstraps(
            n_draws, coef_bootstraps, cov_bootstraps)

        return coef_draws