程序猿成长史(一):初探自生成数据,最小二乘法线性拟合及非线性多项式拟合

近来刚好在实验室里,学习的过程中刚好碰到了人工智能最基础的方面,线性拟合。同时也是接到实验室里一个大佬的任务,生成所需线性拟合的数据集。然后也就顺手整理写下了这篇文章。主要内容包括:

  • 数据生成
  • 基于最小二乘法的线性拟合
  • 基于梯度下降+最小二乘法的非线性多项式拟合

然后这次,博主也是很认真的用了Python来实现整一个过程。(认真脸)
话不多说,直接开始吧。

  1. 数据生成
    首先,数据的生成是自己写的函数,思路大概如下:首先是选定x和y的大致范围和生成的数据的个数(在我的代码中,x和y都是一百个数据,x的范围是从0到100,y的范围是从0到200),然后就是选取合适的函数,再逐一生成即可。Talk is cheap , let’s see the code.
import Random
Datax = []
Datay = []	
noise = 10 				#噪声
FilePath="/XXXXX/XXXX/XXXXX/文件储存/4.txt" #自己想保存的文件路径
while(len(Datay)<100):          #生成一百个数据点
    x=random.uniform(0,100)
    y=x*2+6+random.uniform(-noise, noise) 		#生成y=2x+6的函数的点
    if(y<=200):				#选择其中y小于200的点
        Datay.append(y)
        Datax.append(x)

File = open(FilePath,mode='w')			#打开文件并以txt的格式写入
for i in range(100):
    File.write(str(Datax[i])+'\t'+str(Datay[i])+'\n')
    
# 用Excel表保存的写法
# workbook = xlwt.Workbook()
# sheet1 = workbook.add_sheet('DataSet',cell_overwrite_ok=True)
# for i in range(100):
#     sheet1.write(i,0,str(Datax[i]))
#     sheet1.write(i,1,str(Datay[i]))
# workbook.save(FilePath)

File.close()

保存数据的方法有两种,txt方式和excel方式,大家可以自行选择。 但是要注意的是,在生成的数据点的时候,要加上适当的噪声,以便更好的反应你的算法的泛化性。同时有心的可以生成一些随机的噪声点,然后可以通过一些方法去除这些对算法正确性干扰较大的点。我其实是做了五个噪声点的但是发现其实并没有特别好的去处纯噪点的方法,所以又舍去了 然后这个算法生成的数据都是乱序的,大家可以事先sorted排序点,再用matplotlib去检验数据点的准确性。

  1. 基于最小二乘法的线性拟合
    线性拟合这部分,我是参考着西瓜书《机器学习》里面的公式,自己写的方法,然后将预测的直线也可视化了。用的最小二乘法其实也就是求出预测函数和实际的y的值的平方差,求他的最小值,也就是y=wt+b中对w和b求偏导,并令它们偏导为零的情况。算法证明等大家可以自行百度一下,毕竟我也只是初次探索,未必有别家大牛理解那么深刻,先上程序。
import matplotlib.pyplot as plt
import numpy as np

def ReadFile(Filepath,length):			#从txt.文件中读取dataset
    file = open(Filepath,'r')
    datax = []
    datay = []

    for i in range(length):
        line = file.readline().replace('\n','')
        tem=line.split('\t')
        x = float(tem[0])
        y = float(tem[1])
        datax.append(x)
        datay.append(y)
    return datax,datay

def Cal(datax,datay,length):
    avgx = np.mean(datax)
    tem = 0
    squ = 1
    for i in range(length):	#计算wx+b中的w
        tem = tem + datay[i]*(datax[i]-avgx)
        squ = squ + datax[i] * datax[i]

    squ = squ - sum(datax)**2/length
    answ = tem / squ
    b = 0
    for i in range(length):
        b = b + datay[i] - answ * datax[i]
    b = b / length
    return answ,b

def DrawPicture(datax,datay,w,b,c):
	#绘制相应的拟合直线线
    plt.scatter(datax,datay)
    x = [i for i in range(100)]
    y = [i*w+b for i in range(100)]
    plt.plot(x, y, color=c)
    plt.show()

def NumpyRegression(datax,datay):					
	#用了np库中的拟合函数,相比较性能
    w = np.polyfit(datax, datay, 1)
    print(w[0], w[1])
    return w[0], w[1]

if __name__=="__main__":
    Filepath = '/XXXXX/XXXXX/Documents/文件储存/1.txt'
    length = 105
    datax, datay = ReadFile(Filepath, length)
    w, b = Cal(datax, datay, length)
    w1, b1 = NumpyRegression(datax, datay)
    DrawPicture(datax, datay, w, b, 'r') 	#绘制自写预测函数的图
    DrawPicture(datax, datay, w1, b1, 'b')	#绘制np库中函数预测图像

然后代码最后实现的结果大致如下图:

  1. 基于梯度下降+最小二乘法的非线性多项式拟合
    相比较线性拟合,非线性多项式拟合就显得更为的复杂。博主也是从网上看了别人的教程一点点学着写的。网址如下:
    最小二乘法–多项式拟合非线性函数

    代码方面我也借鉴了很多,只有学习率(LearningRate)和训练轮次(iternum)是我后期依据自己的训练集修改的。不得不说梯度下降方法真的是一个靠运气和调参的算法,调参调的好,拟合就好,但是博主最后也没有找到最合适的学习率,故拟合的效果也没有非常非常好。各位可以自己试着调参。

def CalLoss(datax,datay,A):#计算损失函数
    error = 0.0
    for i in range(len(datax)):
        y1 = 0.0
        for k in range(len(A)):
            y1 += A[k]*(datay[i]**k)
        error += (datay[i]-y1) ** 2
        #计算从0到最大值n的所有误差
    return error

def FittingCurveGradient(datax, datay, order, iternum= 41100, learnrate=0.0000000000001):#order所求多项式最高次数
    A = [0.0]*(order+1)     #初始化order个0
    tem = A
    for r in range(iternum+1):
        for k in range(len(A)):
            gradient = 0.0
            for i in range(len(datax)):#100个
                y1 = 0.0
                for j in range(len(A)):#计算每一轮ak总的梯度值
                    y1 += A[j] * datax[i]**j
                gradient += -2 * (datay[i]-y1) * (datax[i]**k)   #计算A[k]的梯度
            A[k] = A[k] - (learnrate*gradient)  #更新A[k]梯度
            error=CalLoss(datax,datay,A)#用于回溯
            if(r==0):
                errormin=error
            elif(errormin > error):
                errormin=error
                tem[k] = A[k]
            else:
                A[k] = tem[k] #实现回溯
        #检查误差
        if r % 100 == 0:
            error = CalLoss(datax,datay,A)
            print('第{}次迭代,误差下降至:{}'.format(r, error))
    return A

然后自己也偷偷调用了一下np库里的函数去拟合,效果也没有很好。可能是我了解的不够深刻,希望大家和我一起进步。这就是本次项目的基本内容了。期待着未来自己回头看的时候,可以有更多收获嘻嘻。