text: # https://scikit-learn.org/stable/auto_examples/ensemble/plot_bias_variance.html#sphx-glr-auto-examples-ensemble-plot-bias-variance-py
  print(__doc__)
  
  # Author: Gilles Louppe <g.louppe@gmail.com>
  # License: BSD 3 clause
  
  import numpy as np
  import matplotlib.pyplot as plt
  
  from sklearn.ensemble import BaggingRegressor
  from sklearn.tree import DecisionTreeRegressor
  
  # Settings
  n_repeat = 50  # Number of iterations for computing expectations
  n_train = 50  # Size of the training set
  n_test = 1000  # Size of the test set
  n_noise = 0.1  # Standard deviation of the n_noise
  np.random.seed(0)
  
  # Change this for exploring the bias-variance decomposition of other
  # estimators. This should work well for estimators with high variance (e.g.
  # decision trees or KNN), but poorly for estimators with low variance (e.g.
  # linear models)
  estimators = [
      ("Tree", DecisionTreeRegressor()),
      ("Bagging(Tree)", BaggingRegressor(DecisionTreeRegressor())),
  ]
  
  n_estimators = len(estimators)
  
  
  # Generate data
  def f(x):
      x = x.ravel()
  
      return np.exp(-(x ** -2)) + 1.5 * np.exp(-((x - 2) ** 2))
  
  
  def generate(n_samples, n_noise, n_repeat=1):
      X = np.random.rand(n_samples) * 10 - 5
      X = np.sort(X)
  
      if n_repeat == 1:
          y = f(X) + np.random.normal(0.0, n_noise, n_samples)
      else:
          y = np.zeros((n_samples, n_repeat))
  
          for i in range(n_repeat):
              y[:, i] = f(X) + np.random.normal(0.0, n_noise, n_samples)
  
      X = X.reshape((n_samples, 1))
  
      return X, y
  
  
  X_train = []
  y_train = []
  
  for i in range(n_repeat):
      X, y = generate(n_samples=n_train, n_noise=n_noise)
      X_train.append(X)
      y_train.append(y)
  
  X_test, y_test = generate(n_samples=n_test, n_noise=n_noise, n_repeat=n_repeat)
  
  plt.figure(figsize=(10, 8))
  
  # Loop over estimators to compare
  for n, (name, estimator) in enumerate(estimators):
      # Compute predictions
      y_predict = np.zeros((n_test, n_repeat))
  
      for i in range(n_repeat):
          estimator.fit(X_train[i], y_train[i])
          y_predict[:, i] = estimator.predict(X_test)
  
      # Bias^2 + Variance + n_noise decomposition of the MSE
      y_error = np.zeros(n_test)
  
      for i in range(n_repeat):
          for j in range(n_repeat):
              y_error += (y_test[:, j] - y_predict[:, i]) ** 2
  
      y_error /= n_repeat * n_repeat
  
      y_n_noise = np.var(y_test, axis=1)
      y_bias = (f(X_test) - np.mean(y_predict, axis=1)) ** 2
      y_var = np.var(y_predict, axis=1)
  
      print(
          "{0}: {1:.4f} (error) = {2:.4f} (bias^2) "
          " + {3:.4f} (var) + {4:.4f} (n_noise)".format(
              name, np.mean(y_error), np.mean(y_bias), np.mean(y_var), np.mean(y_n_noise)
          )
      )
  
      # Plot figures
      plt.subplot(2, n_estimators, n + 1)
      plt.plot(X_test, f(X_test), "b", label="$f(x)$")
      plt.plot(X_train[0], y_train[0], ".b", label="LS ~ $y = f(x)+noise$")
  
      for i in range(n_repeat):
          if i == 0:
              plt.plot(X_test, y_predict[:, i], "r", label=r"$\^y(x)$")
          else:
              plt.plot(X_test, y_predict[:, i], "r", alpha=0.05)
  
      plt.plot(X_test, np.mean(y_predict, axis=1), "c", label=r"$\mathbb{E}_{LS}\^y(x)$")
  
      plt.xlim([-5, 5])
      plt.title(name)
  
      if n == n_estimators - 1:
          plt.legend(loc=(1.1, 0.5))
  
      plt.subplot(2, n_estimators, n_estimators + n + 1)
      plt.plot(X_test, y_error, "r", label="$error(x)$")
      plt.plot(X_test, y_bias, "b", label="bias^2(x)$")
      plt.plot(X_test, y_var, "g", label="variance(x)$")
      plt.plot(X_test, y_n_noise, "c", label="noise(x)$")
  
      plt.xlim([-5, 5])
      plt.ylim([0, 0.1])
  
      if n == n_estimators - 1:
          plt.legend(loc=(1.1, 0.5))
  
  plt.subplots_adjust(right=0.75)
  plt.show()
  