boot2_slopes = []
boot2_interc = []
n_boots = 100
plt.figure()
for _ in range(n_boots):
 # create a sampling of the residuals with replacement
 boot_resids = np.random.choice(resids, n_points, replace=True)
 y_temp = [y_pred_i + resid_i for y_pred_i, resid_i in zip(y_pred, boot_resids)]

 sample_df = pd.DataFrame({'x': list(x), 'y': y_temp})
 # Fit a linear regression
 ols_model_temp = sm.ols(formula = 'y ~ x', data=sample_df)
 results_temp = ols_model_temp.fit()

 # get coefficients
 boot2_interc.append(results_temp.params[0])
 boot2_slopes.append(results_temp.params[1])

 # plot a greyed out line
 y_pred_temp = ols_model_temp.fit().predict(sample_df['x'])
 plt.plot(sample_df['x'], y_pred_temp, color='grey', alpha=0.2)