<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def _bootstrap_samples_of_smoothing(self, X, y, weights=None, n_bootstraps=1, objective='auto'):
 """Sample the smoothing parameters using simulated response data. For now, the grid of `lam` values is 11 random points in M-dimensional space, where M = the number of lam values, ie len(flatten(gam.lam)) all values are in [1e-3, 1e3] """

        mu = self.predict_mu(X)  # Wood pg. 198 step 1
        coef_bootstraps = [self.coef_]
        cov_bootstraps = [
            load_diagonal(self.statistics_['cov'])]

        for _ in range(n_bootstraps - 1):  # Wood pg. 198 step 2
            # generate response data from fitted model (Wood pg. 198 step 3)
            y_bootstrap = self.distribution.sample(mu)

            # fit smoothing parameters on the bootstrap data
            # (Wood pg. 198 step 4)
            # TODO: Either enable randomized searches over hyperparameters
            # (like in sklearn's RandomizedSearchCV), or draw enough samples of
            # `lam` so that each of these bootstrap samples get different
            # values of `lam`. Right now, each bootstrap sample uses the exact
            # same grid of values for `lam`, so it is not worth setting
            # `n_bootstraps > 1`.
            gam = deepcopy(self)
            gam.set_params(self.get_params())

            # create a random search of 11 points in lam space
            # with all values in [1e-3, 1e3]
            lam_grid = np.random.randn(11, len(flatten(self.lam))) * 6 - 3
            lam_grid = np.exp(lam_grid)
            gam.gridsearch(X, y_bootstrap, weights=weights, lam=lam_grid,
                           objective=objective)
            lam = gam.lam

            # fit coefficients on the original data given the smoothing params
            # (Wood pg. 199 step 5)
            gam = deepcopy(self)
            gam.set_params(self.get_params())
            gam.lam = lam
            gam.fit(X, y, weights=weights)

            coef_bootstraps.append(gam.coef_)

            cov = load_diagonal(gam.statistics_['cov'])

            cov_bootstraps.append(cov)
        return coef_bootstraps, cov_bootstraps