python中怎么推导线性回归模型

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import MinMaxScaler
import matplotlib.pyplot as plt

import seaborn as sns
sns.set_context('notebook')
sns.set_style('white')

model_data = pd.read_csv('model_data.csv',engine='python')
model_data.head()

feat_cols = model_data.columns.tolist()[1:]
print(feat_cols)

X =  model_data[feat_cols].values
y = model_data['价格'].values

def Cost_Function(X,y,theta):
    '''
    需要传入的参数为
    X：自变量
    y:应变量
    theta:权重
    使用均方误差（MSE），作为损失函数
    '''
    m = y.size  #求出y的个数（一共多少条数据）
    t = X.dot(theta) #权重和变量点乘，计算出使用当前权重时的预测值
    c = 0 #定义损失值
    c = 1/(2*m) * (np.sum(np.square(t-y))) 
    #预测值与实际值的差值，平方后除以数据的条数，计算出均方误差。最后乘以1/2（无实际意义，方便以后计算）
    return c #返回损失值

def GradientDescent(X, y, feat_cols, alpha=0.3, num_iters=10000):
    '''
    需传入参数为
    X：自变量
    y:应变量
    feat_cols：变量列表
    alpha：学习率，默认0.3
    num_iters：迭代次数，默认10000次
    使用梯度下降法迭代权重
    '''
    scaler = MinMaxScaler() #最大最小值归一化自变量
    X = scaler.fit_transfORM(X) #归一化
    m = y.size #求出y的个数（一共多少条数据）
    J_history = np.zeros(num_iters) #创建容纳每次迭代后损失值得矩阵，初始值为0
    theta = np.zeros(len(feat_cols)+1) #设置默认权重，0
    for iter in np.arange(num_iters): #根据迭代次数，开始迭代
        t = X.dot(theta) #权重和变量点乘，计算出使用当前权重时的预测值
        theta = theta - alpha*(1/m)*(X.T.dot(t-y))
        #对代价函数求导，算出下降最快的方向，乘以学习率（下降的速度），再用原来的权重相减，得到新的权重
        J_history[iter] = Cost_Function(X, y, theta) #求出新的权重时的损失值，存入矩阵
    return(theta, J_history) #返回最终的权重和历次迭代的损失值

def plot_Cost(GD_result):
    '''
    绘制权重变化情况
    需传入参数为
    GD_result：梯度下降法结果
    '''
    theta , Cost = GD_result #得到权重和损失值
    print('theta: ',theta.ravel()) #打印权重

    plt.plot(Cost) #绘制损失值变化情况
    plt.title('COST change')
    plt.ylabel('Cost')
    plt.xlabel('Iterations')
    plt.grid()
    plt.show()

def lr_function(X,y,feat_cols):
    '''
    需要输入的变量为
    X：自变量
    y:应变量
    feat_cols：变量列表
    '''
    def score(y_p,y):
        '''
        y_p:预测值
        y:真实值
        dimension：样本数量
        计算R^和调整R^    
        '''
        aa=y_p.copy(); bb=y.copy()
        if len(aa)!=len(bb):
            print('not same length')
            return np.nan

        cc=aa-bb
        wcpfh=sum(cc**2) #误差平方和

        # RR means R_Square
        RR=1-sum((bb-aa)**2)/sum((bb-np.mean(bb))**2)
        return RR#返回R^
    X = np.c_[np.ones(X.shape[0]),X]
    GD_result = GradientDescent(X, y, feat_cols)
    plot_Cost(GD_result)
    y_p = np.dot(X,GD_result[0])
    RR = score(y_p,y)
    return RR,y_p,GD_result[0]  #返回R^,预测值

model_result = lr_function(X,y,feat_cols)
print('R^2为:{}'.format(round(model_result[0],4)))

scaler = MinMaxScaler()
X = scaler.fit_transform(X)

lr = LinearRegression()
lr.fit(X,y)
R2 = lr.score(X,y)
print('R^2为:{}'.format(round(R2,4)))