()
| 110 | |
| 111 | |
| 112 | def plot_bayes(): |
| 113 | np.random.seed(12345) |
| 114 | n_in = 1 |
| 115 | n_out = 1 |
| 116 | n_ex = 20 |
| 117 | std = 15 |
| 118 | intercept = 10 |
| 119 | X_train, y_train, X_test, y_test, coefs = random_regression_problem( |
| 120 | n_ex, n_in, n_out, intercept=intercept, std=std, seed=0 |
| 121 | ) |
| 122 | |
| 123 | # add some outliers |
| 124 | x1, x2 = X_train[0] + 0.5, X_train[6] - 0.3 |
| 125 | y1 = np.dot(x1, coefs) + intercept + 25 |
| 126 | y2 = np.dot(x2, coefs) + intercept - 31 |
| 127 | X_train = np.vstack([X_train, np.array([x1, x2])]) |
| 128 | y_train = np.hstack([y_train, [y1[0], y2[0]]]) |
| 129 | |
| 130 | LR = LinearRegression(fit_intercept=True) |
| 131 | LR.fit(X_train, y_train) |
| 132 | y_pred = LR.predict(X_test) |
| 133 | loss = np.mean((y_test - y_pred) ** 2) |
| 134 | |
| 135 | ridge = RidgeRegression(alpha=1, fit_intercept=True) |
| 136 | ridge.fit(X_train, y_train) |
| 137 | y_pred = ridge.predict(X_test) |
| 138 | loss_ridge = np.mean((y_test - y_pred) ** 2) |
| 139 | |
| 140 | LR_var = BayesianLinearRegressionKnownVariance( |
| 141 | mu=np.c_[intercept, coefs][0], sigma=np.sqrt(std), V=None, fit_intercept=True, |
| 142 | ) |
| 143 | LR_var.fit(X_train, y_train) |
| 144 | y_pred_var = LR_var.predict(X_test) |
| 145 | loss_var = np.mean((y_test - y_pred_var) ** 2) |
| 146 | |
| 147 | LR_novar = BayesianLinearRegressionUnknownVariance( |
| 148 | alpha=1, beta=2, mu=np.c_[intercept, coefs][0], V=None, fit_intercept=True |
| 149 | ) |
| 150 | LR_novar.fit(X_train, y_train) |
| 151 | y_pred_novar = LR_novar.predict(X_test) |
| 152 | loss_novar = np.mean((y_test - y_pred_novar) ** 2) |
| 153 | |
| 154 | xmin = min(X_test) - 0.1 * (max(X_test) - min(X_test)) |
| 155 | xmax = max(X_test) + 0.1 * (max(X_test) - min(X_test)) |
| 156 | X_plot = np.linspace(xmin, xmax, 100) |
| 157 | y_plot = LR.predict(X_plot) |
| 158 | y_plot_ridge = ridge.predict(X_plot) |
| 159 | y_plot_var = LR_var.predict(X_plot) |
| 160 | y_plot_novar = LR_novar.predict(X_plot) |
| 161 | |
| 162 | y_true = [np.dot(x, coefs) + intercept for x in X_plot] |
| 163 | fig, axes = plt.subplots(1, 4) |
| 164 | |
| 165 | axes = axes.flatten() |
| 166 | axes[0].scatter(X_test, y_test) |
| 167 | axes[0].plot(X_plot, y_plot, label="MLE") |
| 168 | axes[0].plot(X_plot, y_true, label="True fn") |
| 169 | axes[0].set_title("Linear Regression\nMLE Test MSE: {:.2f}".format(loss)) |
nothing calls this directly
no test coverage detected