线性回归 算法 实例

1.000000 0.067732 3.176513
1.000000 0.427810 3.816464
1.000000 0.995731 4.550095
1.000000 0.738336 4.256571
1.000000 0.981083 4.560815
1.000000 0.526171 3.929515
1.000000 0.378887 3.526170
1.000000 0.033859 3.156393
1.000000 0.132791 3.110301
1.000000 0.138306 3.149813
1.000000 0.247809 3.476346
1.000000 0.648270 4.119688
1.000000 0.731209 4.282233
1.000000 0.236833 3.486582
1.000000 0.969788 4.655492
1.000000 0.607492 3.965162
1.000000 0.358622 3.514900
1.000000 0.147846 3.125947
1.000000 0.637820 4.094115
1.000000 0.230372 3.476039

#coding:utf-8
from numpy import *

#获取样本数据:
#前面是特征数据，最后一列是数据的输出
def loadDataSet(fileName):      
    #特征维数
    numFeat = len(open(fileName).readline().split('\t')) - 1 
    #特征值
    dataMat = [];
    #真实目标值(输出值) 
    labelMat = []
    #打开文件
    fr = open(fileName)
    #按行读取
    for line in fr.readlines():
        #构建列表，存储.txt文件每一行数据
        lineArr =[] #list类型
        #分割数据
        curLine = line.strip().split('\t')
        #遍历每一维度
        for i in range(numFeat):
            #该行的每一个维度数值存入到lineArr中
            lineArr.append(float(curLine[i]))
        #存所有行的数据的特征
        dataMat.append(lineArr)
        #存目标值(行向量)
        labelMat.append(float(curLine[-1])) 
    return dataMat,labelMat

#最小二乘法（标准回归函数）：求拟合直线的参数w=(X.T*X).I*X.T*y
#样本特征数据xArr
#样本的目标值yArr
def standRegres(xArr,yArr):
    #转成矩阵类型
    xMat = mat(xArr) 
    #行向量转置为列向量形式
    yMat = mat(yArr).T
    #
    xTx = xMat.T*xMat 
    #对xTx求逆，需要满足其行列式不为0
    if linalg.det(xTx) == 0.0:
        print "This matrix is singular, cannot do inverse"
        return
    #拟合直线的参数
    ws = xTx.I * (xMat.T*yMat)
    return ws

#***批处理梯度下降法
#输入为训练数据与其输出
def gradDescent(dataMatIn, classLabels):
    #转为NumPy可识别的矩阵
    dataMatrix = mat(dataMatIn)             
    #为了便于计算，classLabels为行向量转为列向量
    labelMat = mat(classLabels).transpose() 
    #获取输入数据的条数m，特征数n
    m,n = shape(dataMatrix)    
    #设定迭代的步长alpha
    alpha = 0.001
    #设置循环次数500次，即训练次数，人为给定
    maxCycles = 500
    #权值初始化为1，后面根据样本数据调整
    #训练结束得到最优权值
    #weights为n行，1维。为列向量。
    weights = ones((n,1))
    #循环maxCycles次，
    #每次根据模型输出结果与真实值的误差，调整权值。
    for k in range(maxCycles):    
        #dataMatrix*weights矩阵的乘法。
        #事实上包含600次的乘积   
        #h为模型给出的一个预测值       
        h = dataMatrix*weights  
        #计算误差，每条记录真实值与预测值之差     
        error = h-labelMat
        #权值调整(未知参数调整)        
        weights = weights - alpha * dataMatrix.transpose()* error
    #循环次数结束，返回回归系数  
    return weights

#***随机梯度下降法
#输入：训练数据矩阵、输出、迭代次数
#输出：权值w
def stocGradDescent0(dataMatrix, classLabels, numIter=100):
    #list类型转换array类型
    dataMatrix=array(dataMatrix)
    #获取训练数据条数m，维度n
    m,n = shape(dataMatrix)
    #步长
    alpha = 0.01
    #权值初始化为1，后面根据样本数据调整
    weights = ones(n)   
    #循环numIter次
    for j in range(numIter):
        #遍历每一条数据
        for i in range(m):

        #h为当前样本的预测值，批处理梯度上升算法的h为所有样本的模型输出
        #此处h为一个值，即一次只有一个样本更新
        #dataMatrix[i]*weights也为当前样本行乘以权值weights
            h = sum(dataMatrix[i]*weights)
    #误差，此处error为一个值
            error = h-classLabels[i]
    #只选择当前样本进行权值更新
            #weights = weights - dot(alpha * error,dataMatrix[i])#list类型运算
            weights = weights - alpha * error*dataMatrix[i]#array类型直接数乘
    #返回权值   
    return mat(weights).transpose()


#图形化显示标准线性回归结果，包括数据集及它的最佳拟合直线    
def standplot(xarr,yarr,w):
    import matplotlib.pyplot as plt
    xmat=mat(xarr)
    ymat=mat(yarr)
    #画点
    fig=plt.figure()
    ax=fig.add_subplot(111)
    ax.scatter(xmat[:,1].flatten().A[0],ymat.T[:,0].flatten().A[0])
    #画线，为了保证直线上的点是按顺序排列，需按升序排列
    xCopy = xmat.copy()
    #排序
    xCopy.sort(0)  #排序
    yHat = xCopy*w #预测值
    ax.plot(xCopy[:,1],yHat)
    plt.show() 


if __name__ == "__main__":    
    dataMatIn, classLabels = loadDataSet('ex0.txt')
    #最小二乘解析式法
    w_LS=standRegres(dataMatIn,classLabels)
    print'w_LS='
    print w_LS
    standplot(dataMatIn,classLabels,w_LS)

    #批处理梯度下降法
    w_gradDescent=gradDescent(dataMatIn, classLabels)
    print'w_gradDescent='
    print w_gradDescent
    #standplot(dataMatIn,classLabels,w_gradDescent)

    #随机梯度下降法
    w_stocGradDescent0=stocGradDescent0(dataMatIn, classLabels)
    print'w_stocGradDescent0='
    print w_stocGradDescent0
    #standplot(dataMatIn,classLabels,w_stocGradDescent0)