from astroML.linear_model import LinearRegression
linear = LinearRegression()
linear.fit(z_sample[:,None], mu_sample, dmu)
mu_fit = linear.predict(z[:, None])

#------------------------------------------------------------
# Plot the results
#fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111)
ax.plot(z, mu_fit_sk, '-k')
ax.plot(z, mu_fit, '-k', color='red')
ax.plot(z, mu_true, '--', c='gray')
ax.errorbar(z_sample, mu_sample, dmu, fmt='.k', ecolor='gray', lw=1)
ax.set_xlim(0.01, 1.8)
ax.set_ylim(36.01, 48)
ax.set_ylabel(r'$\mu$')
ax.set_xlabel(r'$z$')
plt.show()