使用Python手撸了一个线性回归模型
## 线性回归模型 import numpy as np import pandas as pd import matplotlib.pyplot as plt def getData(): r""" 使用numpy构造满足条件的随机数() :return: """ dataX = np.random.randint(1,100,size=(1,100)) dataY = np.random.randint(1,100,size=(1,100)) data = np.concatenate((dataX,dataY), axis = 0) G1tfa = data > 50 G1tfi = np.logical_and(G1tfa[0,:], G1tfa[1,:]) G1index = np.nonzero(G1tfi)[0] G1 = data[:, G1index] G2tfa = data < 50 G2tfi = np.logical_and(G2tfa[0,:], G2tfa[1,:]) G2index = np.nonzero(G2tfi)[0] G2 = data[:, G2index] return (G1, G2) def getData2(): r""" 使用numpy生成随机数; 使用pandas构造满足条件的随机数; :return: """ df = pd.DataFrame() df['X'] = np.random.randint(1,100,size=(100)) df['Y'] = np.random.randint(1,100,size=(100)) G1 = df[(df['X']>50) & (df['Y']>50)] G2 = df[(df['X']<50) & (df['Y']<50)] return (G1.values.T, G2.values.T) def getData3(): r""" 生成正太分布的随机数; :return: """ # data = np.random.uniform(0, 1, size = 1000)#随机均匀采样 # data3 = np.random.rand(1000) #随机均匀分布 # # data2 = np.random.normal(0, 1, size = 1000)#正太分布,均值为0,方差为1 # data4 = np.random.randn(1000)#正太分布,均值为0,方差为1 point1 = (0.5, 0.5) point2 = (2.5, 2.5) a = np.random.normal(point1[0], 0.3, size = (1000,1)) b = np.random.normal(point1[1], 0.2, size = (1000,1)) c = np.random.normal(point2[0], 0.2, size = (1000,1)) d = np.random.normal(point2[1], 0.3, size = (1000,1)) G1 = pd.DataFrame(data = np.concatenate((a,b),axis=1), columns=['X','Y']) G2 = pd.DataFrame(data = np.concatenate((c,d),axis=1), columns=['X','Y']) return (G1.values.T, G2.values.T) def getData4(): r""" 生成可以用来线性回归的数据 :return: """ w = 1.2 b = 0.2 x = np.linspace(1,100,100) noise = np.random.randn(100)*10 y = w * x + b + noise return (x,y) def linefit(x, y): r""" 采用梯度下降的方法更新w和b :param x: :param y: :return: """ #初始化迭代参数; w = 0 b = 0 w_last = 0.001 b_last = 0.001 loss_last = 6330 stop = 10000 i = 0 eta = 0.00001#学习率,如果学习率很大,将无法收敛,loss震荡很大; count = 0 #如果100次都没进步,就停止迭代 loss_list = [] #用来保存loss曲线 w_list = [0.001, 0]#用来保存每次迭代后的w b_list = [0.001, 0]#用来保存每次迭代后的b while(i < stop): loss = getLoss(w,b,x,y) w_t = w - eta*(loss-loss_last)/(w-w_last) b_t = b - eta*(loss-loss_last)/(b-b_last) w_list.append(w_t) b_list.append(b_t) w_last = w w = w_t b_last = b b = b_t if abs(loss_last - loss)< 0.01: count += 1 if count >=100: print("w is {:.2f}, b is {:.2f}, error is {}, i is {}".format(w,b,loss,i)) break print("w is {:.2f}, b is {:.2f}, loss is {:.2f},i is {}".format(w,b,loss,i)) loss_list.append(loss) loss_last = loss i += 1 # xb = np.mean(x) # yb = np.mean(y) return w,b,loss_list,w_list,b_list def getLoss(w,b,x,y,MAE = True): if MAE: #L1 mean absolute error (MAE) totalLoss = np.sum(abs((w*x+b)-y)) else: #L2 mean squari error (MSE) totalLoss = np.sum((w*x+b)-y**2) return totalLoss def getErrorSurface(x,y): r""" 获取error surface :param x: :param y: :return: """ minError = 999999 w_f = 0 b_f = 0 size = 200 r = np.zeros((size,size)) for i,b in enumerate(np.linspace(-2.1,2.2,size)): for j,w in enumerate(np.linspace(-2.1,2.2,size)): r[i,j] = getLoss(w,b,x,y) if r[i,j] < minError: minError = r[i,j] w_f = w b_f = b return r,w_f, b_f def getErrorSurface2(x,y): r""" 获取error surface (塞给contour函数) :param x: :param y: :return: """ minError = 999999 w_f = 0 b_f = 0 size = 200 r = np.zeros((size,size)) wa = np.linspace(-2.1,2.2,size) ba = np.linspace(-2.1,2.2,size) z = [] for i,b in enumerate(np.linspace(-2.1,2.2,size)): for j,w in enumerate(np.linspace(-2.1,2.2,size)): r[i,j] = getLoss(w,b,x,y) z.append(r[i,j]) if r[i,j] < minError: minError = r[i,j] w_f = w b_f = b z = np.array(z) za = z.reshape((size,size)) return wa,ba,za if __name__ == "__main__": print("线性回归模型") np.random.seed (10) # 绘制散点图; fig, ax = plt.subplots() # G1, G2 = getData3() # ax.scatter(G1[0,:], G1[1,:]) # ax.scatter(G2[0,:], G2[1,:]) x,y = getData4() w1,b1,lossList,w_list,b_list = linefit(x, y) r,w,b= getErrorSurface(x,y) w2,b2,z2= getErrorSurface2(x,y) ax.scatter(x, y,color = "C0") ax.plot(np.r_[0,100],1.2*np.r_[0,100]+0.2,color = "C1",linewidth=3.0,label='Ground Truth') ax.plot(np.r_[0,100],w*np.r_[0,100]+b,color = "C2",linewidth=3.0,label='fit from error surface') ax.plot(np.r_[0,100],w1*np.r_[0,100]+b1,color = "C3",linewidth=3.0,label='Grident') ax.legend() ax.set_title("GT and scatter ") fig.show() #绘制error surface fig, ax = plt.subplots() ax.contourf(w2,b2,z2) ax.scatter(np.array(w_list), np.array(b_list)) ax.scatter(np.array(1.2), np.array(0.2), color = 'r') ax.annotate('', xy=(1.2, 0.2), xycoords='data', xytext=(1.2+0.4, 0.2+0.4), textcoords='axes fraction', arrowprops=dict(facecolor='white', shrink=0.05), horizontalalignment='right', verticalalignment='top', ) # im = ax.imshow(r) # cb = plt.colorbar(im,) ax.set_xlabel('w') # ax.set_xlim(xmin = 0.2, xmax = 2.2) ax.set_ylabel('b') # ax.set_ylim(ymin = 0, ymax = 0.4) ax.set_title("error surface ") fig.show() # 绘制loss # lossList fig, ax = plt.subplots() ax.plot(lossList) ax.set_xlabel('step') # ax.set_xlim(xmin = 0.2, xmax = 2.2) ax.set_ylabel('loss') ax.set_ylim(ymin = 700, ymax = 1000) ax.set_title("loss ") fig.show() #绘制直方图; # ax.hist(data,bins=20) #随机均匀采样 # ax.hist(data3,bins=20) #随机均匀分布 # ax.hist(data2,bins=20) #正太分布,均值为0,方差为1 # ax.hist(data4,bins=20) #正太分布,均值为0,方差为1 #固定X轴、Y轴的范围 # ax.set_ylim(ymin = 0, ymax = 130) # ax.set_xlim(xmin = -5, xmax = 5)