python写一个一元线性回归模型，并作图

python写一个一元线性回归模型，并作图。

设方程为：y=a+bx

from scipy import stats
import numpy as np
import pylab

x = np.array([2899.9446,2923.2775,2937.3601,2944.541,2941.007,2952.3383,2932.5058,2908.766,2867.8376,2821.4957,2777.5559,2768.6795,2794.5523,2774.7532,2814.9944,2797.2596,2808.9146,2815.7976,2823.8238,2883.096,2880.0012,2880.3307,2883.4358,2897.4253,2863.5673,2902.1932,2893.7564,2890.9192,2886.2365,2924.1063,2930.1538,2957.4116,2985.8648,2999.6013,3024.7388,3021.2024,3008.8118,3031.2351,3030.7544,2978.1178,2985.6586,2999.2789,3006.4467,2977.0767,2985.3406,2955.4325,2929.0875,2932.167,2905.1892,2913.5704])
y = np.array([21.65,21.59,21.97,22.1,22.06,21.98,21.99,21.99,21.7,21.72,21.7,21.69,21.35,21.17,21.6,21.5,21.45,21.8,22.2,23.45,24.09,24.65,25.44,26.34,25.85,26.64,26.12,26.39,26.1,27.98,30.55,30.66,32.38,32,30.84,30.9,29.99,29.76,29.87,29.8,29.59,29.88,29.79,30.53,30.58,29.53,29.75,29.78,29.4,30.06])

slope, intercept, r_value, p_value, slope_std_error = stats.linregress(x, y)

predict_y = intercept + slope * x
pred_error = y - predict_y
degrees_of_freedom = len(x) - 2
residual_std_error = np.sqrt(np.sum(pred_error**2) / degrees_of_freedom)
print('相关系数Rxy:',r_value)
print('截距--参数a：',intercept)
print('斜率--参数b：',slope)
# Plotting
pylab.plot(x, y, 'o')
pylab.plot(x, predict_y, 'k-')
pylab.show()