x0 = np.arange(-5, 5)
y0 = alpha_in + beta_in[0] * x0

from sklearn.linear_model import LinearRegression as LinearRegression_sk
from astroML.linear_model import LinearRegression
from astroML.linear_model import LinearRegressionwithErrors

linreg_sk = LinearRegression_sk()
linreg_sk.fit(xi[0][:, None], yi)
linreg_sk_y_fit = linreg_sk.predict(x0[:, None])

linreg = LinearRegression()
linreg.fit(xi[0][:, None], yi, yi_error)
linreg_y_fit = linreg.predict(x0[:, None])

linreg_xy_err = LinearRegressionwithErrors()
linreg_xy_err.fit(xi, yi, yi_error, xi_error)

linreg_xy_err_y_fit = linreg_xy_err.coef_[0] + linreg_xy_err.coef_[1] * x0

# Plot the results
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111)
ax.plot(x0, linreg_sk_y_fit, '-k', color='grey', label='sklearn, no errors')
ax.plot(x0, linreg_y_fit, '-k', color='blue', label='astroML y errors only')
ax.plot(x0, linreg_xy_err_y_fit, '-k', color='red', label='astroML x and y errors')
ax.plot(x0, y0, '--', c='black')
ax.errorbar(xi[0], yi, yi_error, xi_error[0], fmt='.k', ecolor='gray', lw=1)

ax.set_ylabel(r'$y$')
ax.set_xlabel(r'$x$')
ax.legend()

plt.show()