Simulate from the posterior of the coefficients and smoothing params.

        Samples are drawn from the posterior of the coefficients and smoothing
        parameters given the response in an approximate way. The GAM must
        already be fitted before calling this method; if the model has not
        been fitted, then an exception is raised. Moreover, it is recommended
        that the model and its hyperparameters be chosen with `gridsearch`
        (with the parameter `keep_best=True`) before calling `sample`, so that
        the result of that gridsearch can be used to generate useful response
        data and so that the model's coefficients (and their covariance matrix)
        can be used as the first bootstrap sample.

        These samples are drawn as follows. Details are in the reference below.

        1. ``n_bootstraps`` many "bootstrap samples" of the response (``y``) are
        simulated by drawing random samples from the model's distribution
        evaluated at the expected values (``mu``) for each sample in ``X``.

        2. A copy of the model is fitted to each of those bootstrap samples of
        the response. The result is an approximation of the distribution over
        the smoothing parameter ``lam`` given the response data ``y``.

        3. Samples of the coefficients are simulated from a multivariate normal
        using the bootstrap samples of the coefficients and their covariance
        matrices.

        Notes
        -----
        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)
              empirical input data

        y : array of shape (n_samples,)
              empirical response vector

        quantity : {'y', 'coef', 'mu'}, default: 'y'
            What quantity to return pseudorandom samples of.
            If `sample_at_X` is not None and `quantity` is either `'y'` or
            `'mu'`, then samples are drawn at the values of `X` specified in
            `sample_at_X`.

        sample_at_X : array of shape (n_samples_to_simulate, m_features) or
        None, optional
            Input data at which to draw new samples.

            Only applies for `quantity` equal to `'y'` or to `'mu`'.
            If `None`, then `sample_at_X` is replaced by `X`.

        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=5)
            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 parameter is used, and the distribution over the
            smoothing parameters is not estimated using bootstrap sampling.

        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
        -------
        draws : 2D array of length n_draws
            Simulations of the given `quantity` using samples from the
            posterior distribution of the coefficients and smoothing parameter
            given the response data. Each row is a pseudorandom sample.

            If `quantity == 'coef'`, then the number of columns of `draws` is
            the number of coefficients (`len(self.coef_)`).

            Otherwise, the number of columns of `draws` is the number of
            rows of `sample_at_X` if `sample_at_X` is not `None` or else
            the number of rows of `X`.

        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(self, X, y, quantity='y', sample_at_X=None,
               weights=None, n_draws=100, n_bootstraps=5, objective='auto'):
        """Simulate from the posterior of the coefficients and smoothing params.

        Samples are drawn from the posterior of the coefficients and smoothing
        parameters given the response in an approximate way. The GAM must
        already be fitted before calling this method; if the model has not
        been fitted, then an exception is raised. Moreover, it is recommended
        that the model and its hyperparameters be chosen with `gridsearch`
        (with the parameter `keep_best=True`) before calling `sample`, so that
        the result of that gridsearch can be used to generate useful response
        data and so that the model's coefficients (and their covariance matrix)
        can be used as the first bootstrap sample.

        These samples are drawn as follows. Details are in the reference below.

        1. ``n_bootstraps`` many "bootstrap samples" of the response (``y``) are
        simulated by drawing random samples from the model's distribution
        evaluated at the expected values (``mu``) for each sample in ``X``.

        2. A copy of the model is fitted to each of those bootstrap samples of
        the response. The result is an approximation of the distribution over
        the smoothing parameter ``lam`` given the response data ``y``.

        3. Samples of the coefficients are simulated from a multivariate normal
        using the bootstrap samples of the coefficients and their covariance
        matrices.

        Notes
        -----
        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)
              empirical input data

        y : array of shape (n_samples,)
              empirical response vector

        quantity : {'y', 'coef', 'mu'}, default: 'y'
            What quantity to return pseudorandom samples of.
            If `sample_at_X` is not None and `quantity` is either `'y'` or
            `'mu'`, then samples are drawn at the values of `X` specified in
            `sample_at_X`.

        sample_at_X : array of shape (n_samples_to_simulate, m_features) or
        None, optional
            Input data at which to draw new samples.

            Only applies for `quantity` equal to `'y'` or to `'mu`'.
            If `None`, then `sample_at_X` is replaced by `X`.

        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=5)
            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 parameter is used, and the distribution over the
            smoothing parameters is not estimated using bootstrap sampling.

        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
        -------
        draws : 2D array of length n_draws
            Simulations of the given `quantity` using samples from the
            posterior distribution of the coefficients and smoothing parameter
            given the response data. Each row is a pseudorandom sample.

            If `quantity == 'coef'`, then the number of columns of `draws` is
            the number of coefficients (`len(self.coef_)`).

            Otherwise, the number of columns of `draws` is the number of
            rows of `sample_at_X` if `sample_at_X` is not `None` or else
            the number of rows of `X`.

        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 quantity not in {'mu', 'coef', 'y'}:
            raise ValueError("`quantity` must be one of 'mu', 'coef', 'y';"
                             " got {}".format(quantity))

        coef_draws = self._sample_coef(
            X, y, weights=weights, n_draws=n_draws,
            n_bootstraps=n_bootstraps, objective=objective)

        if quantity == 'coef':
            return coef_draws

        if sample_at_X is None:
            sample_at_X = X

        linear_predictor = self._modelmat(sample_at_X).dot(coef_draws.T)
        mu_shape_n_draws_by_n_samples = self.link.mu(
            linear_predictor, self.distribution).T
        if quantity == 'mu':
            return mu_shape_n_draws_by_n_samples
        else:
            return self.distribution.sample(mu_shape_n_draws_by_n_samples)