MLiA笔记_regression

#-*-coding:utf-8-*-

# 8.1 标准回归函数和数据导入函数
from numpy import *

# 用该函数打开一个用tab键分隔的文件文本，默认文件每行的最后一个值是目标值
def loadDataSet(fileName):
    numFeat = len(open(fileName).readline().split('\t')) - 1
    dataMat = []
    labelMat = []
    fr = open(fileName)
    for line in fr.readlines():
        lineArr = []
        curLine = line.strip().split('\t')
        for i in range(numFeat):
            lineArr.append(float(curLine[i]))
        dataMat.append(lineArr)
        labelMat.append(float(curLine[-1]))
    return dataMat, labelMat


# standRegres()用来计算最佳拟合直线
def standRegres(xArr, yArr):
    # 首先读入x和y,并将它们保存到矩阵中；
    xMat = mat(xArr)
    yMat = mat(yArr).T
    # 然后计算xTx，判断它的行列式是否为零
    xTx = xMat.T*xMat
    # numpy提供一个线性代数库linalg，可以直接调用linagle.det()来计算行列式
    if linalg.det(xTx) == 0.0:
        print("this matirx is singular , cannot do inverse")
        return
    # 最后，如果行列式非零，计算并返回w（如果没有检查行列式是否为零就计算矩阵的逆，会出现错误）
    ws = xTx.I*(xMat.T*yMat)
    return ws

# 8.2 局部加权线性回归函数：给定x空间中的任意一点，计算出对应的预测值yHat
# 函数lwlr()，读入数据并创建所需矩阵，之后创建对角加权矩阵weights
def lwlr(testPoint, xArr, yArr, k=1.0):
    xMat = mat(xArr)
    yMat = mat(yArr).T
    m = shape(xMat)[0]
    weights = mat(eye((m)))
    # 算法遍历数据集，计算每个样本点对应的权重值
    for j in range(m):
        # 随着样本点与待遇测点距离的递增，权重将以指数级衰减（输入参数k控制衰减的速度
        diffMat = testPoint - xMat[j,:]
        weights[j,j] = exp(diffMat*diffMat.T/(-2.0*k**2))
    xTx = xMat.T * (weights*xMat)
    if linalg.det(xTx) == 0.0:
        print("this matrix is singular, cannot do inverse")
        return
    ws = xTx.I*(xMat.T*(weights*yMat))
    return testPoint * ws

# 函数lwlrTest()，用于为数据集中每个点调用lwlr()，这有助于求解k的大小
def lwlrTest(testArr,xArr,yArr,k=1.0):
    m = shape(testArr)[0]
    yHat = zeros(m)
    for i in range(m):
        yHat[i] = lwlr(testArr[i],xArr,yArr,k)
    return yHat

# 8.3 岭回归
# 计算回归系数
# 实现了给定λ下的岭回归求解（如果没指定λ，则默认为0.2；lambda是关键字）
def ridgeRegres(xMat, yMat, lam = 0.2):
    # 首先构造矩阵xTx
    # λ设定为0时一样可能产生错误，仍需做一个检查
    xTx = xMat.T*xMat
    denom = xTx + eye(shape(xMat)[1])*lam
    # 最后如果矩阵非奇异，就计算回归系数并返回
    if linalg.det(denom) == 0.0:
        print("this matrix is sigunlar, 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)
    numTestPts = 30
    wMat = zeros((numTestPts,shape(xMat)[1]))
    for i in range(numTestPts):
        ws = ridgeRegres(xMat,yMat,exp(i-10))
        wMat[i,:] = ws.T
    return wMat

# 向前进逐步线性回归
def regularize(xMat):
    inMat = xMat.copy()
    inMeans = mean(inMat,0)
    inVar = var(inMat,0)
    inMat = (inMat - inMeans)/inVar
    return inMat

def rssError(yArr,yHatArr):
    return ((yArr-yHatArr)**2).sum()

# 逐步线性回归算法的实现，它与lasso做法相近但计算简单
# 输入包括：输入数据xArr、预测变量yArr，表示每次迭代需调整的步长参数eps、表示迭代次数的参数numlIt
def stageWise(xArr, yArr, eps = 0.01, numIt = 100):
    #首先将输入的数据转换并存入矩阵中
    xMat = mat(xArr)
    yMat = mat(yArr).T
    # 将特征值按照均值为0，方差为1进行标准化处理
    yMean = mean(yMat,0)
    yMat = yMat - yMean
    xMat = regularize(xMat)
    '''
    xMeans = mean(xMat,0)
    xVar = var(xMat,0)
    xMat = (xMat - xMeans)/xVar
    '''
    m,n = shape(xMat)
    returnMat = zeros((numIt,n))
    # 创建向量ws来保存w的值
    ws = zeros((n,1))
    # 为实现贪心算法建立ws的两个副本
    wsTest = ws.copy()
    wsMax = ws.copy()
    # 优化过程需要迭代numIt次，并且在厄米慈迭代时都打印w向量，用于分析算法执行的过程和效果
    for i in range(numIt):
        print ws.T
        # 该误差初始值为正无穷，经过与所有的误差比较后取最小的误差
        lowestError = inf;
        # 贪心算法在所有特征上运行两次for循环，分别计算增加或减少特征对误差的影响
        for j in range(n):
            for sign in [-1,1]:
                wsTest = ws.copy()
                wsTest[j] += eps*sign
                yTest = xMat*wsTest
                # 使用的是平方误差，通过函数rssError得到
                rssE = rssError(yMat.A,yTest.A)
                if rssE < lowestError:
                    lowestError = rssE
                    wsMax = wsTest
        ws = wsMax.copy()
        returnMat[i,:] = ws.T
    return returnMat

# 8.5 购物信息的获取函数
# 需导入新模块：time.sleep()、json和urllib2
from time import sleep
import json
import urllib2

def serachForSet(retX, retY, setNum, yr,numPce,origPrc):
    # 一开始要休眠10s，为了防止短时间内有过多的API调用
    sleep(10)
    # 拼接查询的URL字符串，添加API的key和带查询的套装信息
    myAPIstr = 'get from code.google.com'
    searchURL = 'https://www.googleapis.com/shopping/search/v1/public/products?key=%s&country=US&q=lego+%d&alt=json' %(myAPIstr,setNum)
    pg = urllib2.urlopen(searchURL)
    # 打开和解析操作通过json.loads()方法实现，完成后将得到一部字典
    retDict = json.loads(pg.red())
    # 从中找出价格和其他信息
    for i in range(len(retDict['items'])):
        try:
            currItem = retDict['items'][i]
            if currItem['product']['condition'] == 'new':
                newFlag = 1
            else: newFlag = 0
            listOfInv = currItem['product']['inventories']
            # 返回的部分结果是一个产品的数组，将在其上进行循环迭代，判断该产品是否新产品并抽取它的价格
            for item in listOfInv:
                # 不完整的套装也会通过检索结果返回，须将这些信息过滤掉
                # （可以统计描述中的关键词或者用贝叶斯方法判断
                sellingPrice = item['orice']
                # 在此仅用一个简单的启发式方法：如果一个套装的价格比原始价格低一半以上，则认为该套装不完整
                if sellingPrice > origPrc*0.5:
                    print("%d\t%d\t%d\t%f\t%f" % (yr, numPce, newFlag, origPrc, sellingPrice) )
                    # 解析成功后的套装将在屏幕上显示出来并保存在list对象retX和retY中
                    retX.append([yr, numPce, newFlag, origPrc])
                    retY.append((sellingPrice))
        except: print("problem with item %d" %i)

# 负责多次调用searchForSer()
def setDataCollect(retX, retY):
    serachForSet(retX, retY, 8288, 2006, 800, 49.99)
    serachForSet(retX, retY, 10030, 2002, 3096, 269.99)
    serachForSet(retX, retY, 10179, 2007, 5195, 499.99)
    serachForSet(retX, retY, 10181, 2007, 3428, 199.99)
    serachForSet(retX, retY, 10189, 2008, 5922, 299.99)
    serachForSet(retX, retY, 10196, 2009, 3263, 249.99)


# 交叉验证测试岭回归
# crisValidation()函数有三个参数，
# xArr、yArr存有数据集中的X和Y值得list对象，默认lgX、lgY具有相同的长度
# numVal是交叉验证的次数，如果没有指定，取默认值10
def crosValidation(xArr, yArr, numVal = 10):
    # 首先计算数据点的个数m
    m = len(yArr)
    indexList = range(m)
    errorMat = zeros((numVal,30))
    for i in range(numVal):
        # 创建训练集合训练集容器
        trainX = []
        trainY = []
        testX = []
        testY = []
        # 创建一个list并使用random.shuff()函数对其中的元素进行混洗，实现训练集或测试集数据点的随机选取
        random.shuffle(indexList)
        for j in range(m):
            # 将数据集的90%分割成训练集，其余10%为测试集，并将二者分别放入对应容器中
            if j < m*0.9:
                trainX.append(xArr[indexList[j]])
                trainY.append(yArr[indexList[j]])
            else:
                testX.append(xArr[indexList[j]])
                testY.append(yArr[indexList[j]])
    # 对数据混洗后，建立一个新的矩阵wMat来保存岭回归中的所有回归系数
    wMat = ridgeTest(trainX, trainY)
    # 8.4.1中ridgeTest()试用了30个不同的λ值创建了30组不同的回归系数。接下来也要在上述测试集上用30组回归系数来循环测试回归效果
    for k in range(30):
        # 岭回归需要使用标准化后额数据，因此测试数据也需要用测试集相同的参数来执行标准化
        matTestX = mat(testX)
        matTrainX = mat(trainX)
        meanTrain = mean(matTrainX, 0)
        varTrain = var(matTrainX, 0)
        matTestX = (matTestX - meanTrain)/varTrain
        yEst = matTestX * mat(wMat[k,:]).T + mean(trainY)
        # 用rssError()计算误差，并将结果保存在errorMat中
        errorMat[i,k] = rssError(yEst.T.A, arange(testY))
    # 在所有交叉验证完成后，errorMat保存了eidgeTest()里每个λ对应的多个误差值
    meanErrors = mean(errorMat, 0 )
    # 为了将得出的回归系数与standRegres()做对比，需要计算这些误差估计值的均值
    minMean = float(min(meanErrors))
    bestWeights = wMat[nonzero(meanErrors == minMean)]
    # 岭回归使用了数据标准化，而standRegres()则没有，因此为了上述比较须将数据还原
    xMat = mat(xArr)
    yMat = mat(yArr).T
    meanX = mean(xMat, 0)
    varX = var(xMat, 0)
    unReg = bestWeights/varX
    print("the best model from ridge regression is :\n", unReg)
    print("with constant term: "-1*sum(multiply(meanX,unReg)) + mean(yMat))