NumPy实现线性回归

本文主要是介绍NumPy实现线性回归，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

1 单变量线性回归

1.1 sklearn实现(最小二乘法)

import osimport pandas as pd
import matplotlib.pyplot as plt
import syscurrent_dir=os.getcwd()
path=current_dir+'\\'+"Salary Data.csv"def plot_data(path):table=pd.read_csv(path)experience = table["Experience Years"]salary = table["Salary"]plt.figure(figsize=(8,6))plt.scatter(experience,salary,color="blue",label="Data points")plt.title("experience vs year")plt.xlabel("Experience (Years)")plt.ylabel("Salary")plt.grid(True)plt.legend()plt.show()
plot_data(path)table=pd.read_csv(path)
y=table['Salary']
x=table[ ['Experience Years'] ]  # x.shape=(40,1)
z=table['Experience Years']    # z.shape=(40,)from sklearn.model_selection import train_test_split
x_train,x_test,y_train,y_test = train_test_split(x,y,train_size=0.7,random_state=2529)
# (28, 1) (28,) (12, 1) (12,)from sklearn.linear_model import LinearRegression
model = LinearRegression()
model.fit(x_train,y_train)print( model.intercept_ )  # 26596.961311068262
print( model.coef_ )       # [9405.61663234]from sklearn.metrics import mean_squared_error, r2_score
y_pred = model.predict(x_test)mse = mean_squared_error(y_test, y_pred)
print( "mse = ", mse )          # 24141421.671440993
r2 = r2_score(y_test, y_pred)
print( "r2 = ", r2 )            # 0.960233432146844y_whole_pred=model.predict(x)
# x.iloc[:,0]可以写成x, 或者x["Experience Years"]
plt.scatter(x.iloc[:,0],y,color="blue",label="Data points")
plt.plot(x,y_whole_pred,color="red",linewidth=2, label='linear regression')plt.xlabel("Experience (Years)")
plt.ylabel("Salary")
plt.legend()
plt.show()

1.2 NumPy实现(梯度下降法)

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import os
from sklearn.model_selection import train_test_split
import sysdef plot_data(path):table=pd.read_csv(path)experience = table["Experience Years"]salary = table["Salary"]plt.figure(figsize=(8,6))plt.scatter(experience,salary,color="blue",label="Data points")plt.title("experience vs year")plt.xlabel("Experience (Years)")plt.ylabel("Salary")plt.grid(True)plt.legend()plt.show()class MyLinearReg:def __init__(self,lr = 0.01, epochs = 1000):self.lr = lrself.epochs = epochsself.w = Noneself.b = Noneself.loss_history = []def fit(self,X,y):m,n = X.shapeself.w = np.zeros(n)self.b = 0for epoch in range(self.epochs):# x(m,n) * w(n,), numpy广播机制矩阵向量乘法y_pred = X @ self.w + self.b  # y_pred(m,)loss = (y_pred - y)           # loss(m,)dcost_dw = (1/m) * (X.T @ loss)dcost_b = (1/m) *  lossdcost_b = np.sum(dcost_b)self.w -= self.lr * dcost_dwself.b -= self.lr * dcost_bsquare_loss = (y_pred-y)**2mean_loss = np.mean(square_loss)self.loss_history.append(mean_loss)if epoch % 100 == 99 :print(f"Epoch {epoch} loss: {mean_loss}")print("Trainning finished.")print("Final parameters:","Slope w=",self.w," Bias b=",self.b)# Final parameters: Slope w= [9853.19132896]  Bias b= 23780.770014707407def predict(self,X):return X @ self.w + self.bdef get_params(self):return self.w, self.b# plot_data(path)
current_dir=os.getcwd()
path=current_dir+'\\'+"Salary Data.csv"
table=pd.read_csv(path)
x = table["Experience Years"].values # x(40,)
y = table["Salary"].values           # y(40,)
#(32,),(8,)(32,)(8,)
x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2, random_state=42)
# (32,) (32,) (8,) (8,)x_train=x_train.reshape(-1,1)
x_test=x_test.reshape(-1,1)
model = MyLinearReg()
model.fit(x_train, y_train)y_pred = model.predict(x_test)from sklearn.metrics import mean_squared_error, r2_scoremse = mean_squared_error(y_test, y_pred)
print( "mse = ", mse )          # mse =  43053815.910611115
r2 = r2_score(y_test, y_pred)
print( "r2 = ", r2 )           # r2 =  0.9165907194371214X=x.reshape(-1,1)
y_whole_pred=model.predict(X)
# x.iloc[:,0]可以写成x, 或者x["Experience Years"]
plt.scatter(x,y,color="blue",label="Data points")
plt.plot(x,y_whole_pred,color="red",linewidth=2, label='linear regression')plt.xlabel("Experience (Years)")
plt.ylabel("Salary")
plt.legend()
plt.show()

Epoch 99 loss: 111815444.20061775
Epoch 199 loss: 81534511.03025383
Epoch 299 loss: 61760636.04682423
Epoch 399 loss: 48848017.74472436
Epoch 499 loss: 40415896.49608463
Epoch 599 loss: 34909602.800390095
Epoch 699 loss: 31313915.621658318
Epoch 799 loss: 28965881.353634194
Epoch 899 loss: 27432581.973080143
Epoch 999 loss: 26431315.92580659
Trainning finished.
Final parameters: Slope w= [9853.19132896]  Bias b= 23780.770014707407
mse =  43053815.910611115
r2 =  0.9165907194371214

2 多变量线性回归

2.1 sklearn实现(最小二乘法)

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import os
import sysdef draw_PairwiseScatter(x,y):num_features = x.shape[1]plt.figure(figsize=(15,10))for i in range(num_features):plt.subplot(3,5,i+1) # 子图的索引从1开始plt.scatter(x[:,i],y,marker='o', color="green", s=15,alpha=0.5)plt.xlabel("Feature {}".format(i+1))plt.ylabel("Label")plt.title("Featurs {} vs Target".format(i+1))plt.tight_layout()plt.show()
def draw_real_pred(x,y,model):y_pred_whole =  model.predict(x)num_features = x.shape[1]plt.figure( figsize=(15,10) )for i in range(num_features):plt.subplot(3,5,i+1)plt.scatter(x[:,i],y,marker='o',color="green", s=15,alpha=0.5)plt.scatter(x[:,i],y_pred_whole,marker="o", color="red", s=15,alpha=0.5)plt.xlabel("Feature {}".format(i+1))plt.ylabel("Label")plt.title("Featurs {} vs Target".format(i+1))plt.tight_layout()plt.show()current_dir = os.getcwd()
path = current_dir + '\\' + "Boston.csv"
house = pd.read_csv(path)y = house['MEDV']                           #  (506,)
X = house.drop(['MEDV'], axis = 1)   #  (506,13)
X=np.array(X)
y=np.array(y)draw_PairwiseScatter(X,y)from sklearn.linear_model import LinearRegression
model = LinearRegression()from sklearn.model_selection import train_test_split
x_train, x_test, y_train, y_test = train_test_split(X, y, train_size = 0.7, random_state = 2529)
#     (354, 13) (152, 13) (354,) (152,)# Ordinary Least Squares 不是梯度下降，不用标准化数据
# theta = (X.T * X)-1 * X.T * y: 最小二乘法
model.fit(x_train,y_train)
print(model.intercept_)
print(model.coef_)y_pred = model.predict(x_test)from sklearn.metrics import  mean_absolute_error, r2_score
print( "mean_absolute_error(y_pred,y_test):", mean_absolute_error(y_pred,y_test) )print ( model.score(x_test,y_test) )
r2 = r2_score(y_test, y_pred)
print(r2)  # 0.6551914852365524draw_real_pred(X,y,model)

2.2 NumPy实现(梯度下降法)

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import os
import sysdef draw_PairwiseScatter(x,y):num_features = x.shape[1]plt.figure(figsize=(15,10))for i in range(num_features):plt.subplot(3,5,i+1)plt.scatter(x[:,i],y,marker='o', color="green", s=15,alpha=0.5)plt.xlabel("Feature {}".format(i+1))plt.ylabel("Label")plt.title("Featurs {} vs Target".format(i+1))plt.tight_layout()plt.show()def draw_real_pred(x,y,model):y_pred_whole =  model.predict(x)num_features = x.shape[1]plt.figure(figsize=(15,10))for i in range(num_features):plt.subplot(3,5,i+1)plt.scatter(x[:,i],y,marker='o',color="green", s=15,alpha=0.5)plt.scatter(x[:,i],y_pred_whole,marker='o', color="red", s=15,alpha=0.5)plt.xlabel("Feature {}".format(i+1))plt.ylabel("Label")plt.title("Featurs {} vs Target".format(i+1))plt.tight_layout()plt.show()class MultipleLinear:def __init__(self,learning_rate=0.01, epochs=1000):self.learning_rate = learning_rateself.epochs = epochsself. theta = Noneself.cost_history = Nonedef fit(self,X,y):X = np.hstack( ( np.ones((X.shape[0],1)), X ) )m,n = X.shapeself.theta = np.zeros(n)self.cost_history = []for epoch in range(self.epochs):y_pred = X @ self.thetagradient = X.T @ (y_pred - y)self.theta -= self.learning_rate * gradient * (1/m)cost = self.compute_cost(X,y)self.cost_history.append(cost)if epoch % 100 == 99:print(f"Epoch {epoch} cost: {cost}")print("Training complete")print ( self.theta )def predict(self,X):m,n = X.shapeX = np.hstack( (np.ones((m,1)), X) )return  X @ self.thetadef compute_cost(self,X,y):m = X.shape[0]y_pred = X @ self.thetasq_errors = (y_pred - y)**2cost = 1 / (2 * m) * np.sum(sq_errors)return costcurrent_dir = os.getcwd()
path = current_dir + '\\' + "Boston.csv"
house = pd.read_csv(path)y = house['MEDV']                           #  (506,)
X = house.drop(['MEDV'], axis = 1)   #  (506,13)
X=np.array(X)
y=np.array(y)draw_PairwiseScatter(X,y)from sklearn.model_selection import train_test_split
x_train, x_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
# (404, 13) (102, 13) (404,) (102,)from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
x_train_scaled = scaler.fit_transform(x_train)  # (404,13)
x_test_scaled = scaler.transform(x_test)        # (102,13)model = MultipleLinear()
model.fit(x_train_scaled, y_train)y_pred = model.predict(x_test_scaled)from sklearn.metrics import  r2_score
r2 = r2_score(y_test,y_pred)
print("r2 = ",r2)       # r2 =  0.6543244875135051draw_real_pred(X,y,model)