初步岭回归

岭回归

1、 机器学习实战

    在机器学习实战中，是这么引入岭回归的：
如果特征比样本还多（n > m）,也就是说输入的数据矩阵X不是满秩矩阵。非满秩矩阵在求逆时会出错（ $(X^TX) ^{-1}$）为了解决这个问题统计学家引入了这么个概念。岭回归最先用来处理特征数多于样本情况，现在也用于在估计中加入偏差，从而得到更好的估计。这里引入λ来限制W之和，通过引入惩罚项，能减少不重要的参数。
    但是公式具体是怎么个意思，书上并没有讲啊……，心情很烦啊！！

2、 交叉验证

k折交叉验证 k-fold cross validation
将全部训练数据S分成k个不相交的子集，假设S中的训练样例个数为m，那么每一个子集有m/k个训练样例，相应的子集为 {s1,s2,…,sk}；
每次从分好的子集中中拿出一个作为测试集，其它k-1个作为训练集；
根据训练训练出模型或者假设函数；
把这个模型放到测试集上，得到分类率；
计算k次求得的分类率的平均值，作为该模型或者假设函数的真实分类率。
这个方法充分利用了所有样本。但计算比较繁琐，需要训练k次，测试k次。

先贴代码吧

#encoding:utf-8
import numpy as py
from numpy import *
#加载数据
def loadDataSet(filename):
    datMat=[]
    labelMat=[]
    numfat = len(open(filename).readline().split('\t')) - 1

    for line in open(filename).readlines():
        linArr =[]
        curline = line.strip().split('\t')
        for i in range(numfat):
            linArr.append(float(curline[i]))
        datMat.append(linArr)
        labelMat.append(float(curline[-1]))
    return datMat,labelMat
def rigeRes(xMat, yMat,lam = 0.2):
    xTx = xMat.T*xMat
    denom = xTx+eye(shape(xMat)[1])*lam
    if linalg.det(denom) == 0:
        print "this matrix is singular,cannot do inverse"
        return
    ws = denom.I*xMat.T*yMat
    return ws
def ridgeTest(xArr,yArr):
    xMat = mat(xArr)
    yMat = mat(yArr).T
    #数据标准化,将特征值变为同等地位，这样不会因为数字差值最大的属性对计算结果影响最大

    yMean = mean(yMat, 0)
    yMat =yMat - yMean

    xMeans = mean(xMat, 0)
    xVar = var(xMat,0)
    xMat = (xMat - xMeans)/xVar
    numTestPts = 30
    wMat = zeros((numTestPts, shape(xMat)[1]))#存储不同lambda对应的回归系数
    #处理30个不同lambda对应岭回归  通过交叉验证选出最好的lambda最为回归系数
    for i in range(numTestPts):
        ws = rigeRes(xMat, yMat, exp(i - 10))
        wMat[i,:] = ws.T
    return wMat
def rssError(xArr,yArr):
    return ((xArr - yArr)**2).sum()
'''
交叉验证测试岭回归

'''  
def crossValidation(xArr, yArr, numVal = 10):#默认交叉验证次数10
    errorMat=zeros((numVal, 30))#每行 代表每次验证不同lambda 对应回归系数产生的误差，
    xMat = mat(xArr)
    m,n = shape(xMat)
    indexList = range(m)#存储样本个数的下标
    for i in range(numVal):#交叉验证的次数
        #数据分为训练集和测试集容器
        testX=[];trainX=[]#
        testY=[];trainY=[]
        random.shuffle(indexList)  #对其中的元素进行混洗，因此可以进行随机提取
        for j in range(m):
            if j < m * 0.9:#90%为训练集
                trainX.append(xArr[indexList[j]])
                trainY.append(yArr[indexList[j]])
            else:
                testX.append(xArr[indexList[j]])
                testY.append(yArr[indexList[j]])
        ws = ridgeTest(trainX, trainY)#会产生30行n列(n个特征值)de回归系数 具体上面岭回归函数那块
        for k in range(30):#30代表前面30个lambda
            #岭回归需要使用标准化后的数据，因此测试数据也需要用与测试集相同的参数来执行标准化
            matTestX = mat(testX)
            matTrainX = mat(trainX)
            meanTrain = mean(matTrainX, 0)
            varTrain = var(matTrainX, 0)
            matTestX = (matTestX - meanTrain)/varTrain#标准化
            yEst = matTestX*mat(ws[k,:]).T+mean(trainY)#计算误差

            errorMat[i,k] = rssError(yEst.T.A, testY)#第i次验证，第k个lambda对应的误差

    meanErrors = mean(errorMat, 0)#求出不同lambda对应的平均误差值
    minMean = float(min(meanErrors))#找出最小值
    bestWeights = ws[nonzero(meanErrors == minMean)]#找出最小误差值lambda在ws里面所对应的位置
    #对数据进行还原
    xMat = mat(xArr)
    yMat =mat(yArr).T
    meanX = mean(xMat, 0)
    varX = var(xMat, 0)
    unReg = bestWeights/varX
    print "the model from Regression is:\n",unReg#还原后的回归系数
    print "with constant term:"
    print -1*sum(multiply(meanX, unReg)) + mean(yMat)

def show(ridgeWeights):

    import matplotlib.pyplot as plt
    fig = plt.figure()
    ax = fig.add_subplot(111)
#     ax.plot(ridgeWeights)#接收一个矩阵，横坐标就是第几列，纵坐标就是矩阵第一列里面的所有数
    #比如说我要化w0 w1
#     ax.plot(ridgeWeights[:,0]) #所有行第一列的代表w0随lambda值变化而变化
#     ax.plot(ridgeWeights[:,1])#所有行第二列的代表w1随lambda值变化而变化
    plt.show()
abx,aby = loadDataSet('abalone.txt')
ridgeWeights = ridgeTest(abx, aby)
print ridgeWeights
show(ridgeWeights)
crossValidation(abx, aby)