<![CDATA[Code Producers]]>https://codeproducers.com/https://codeproducers.com/favicon.pngCode Producershttps://codeproducers.com/Ghost 3.2Thu, 22 Oct 2026 02:38:39 GMT60<![CDATA[《机器学习实战》笔记 (第八章和第九章) (附Python3版代码)(Machine Learning in Action)]]>https://codeproducers.com/20180402/5e0456874b43c3452b5136deMon, 02 Apr 2018 09:02:01 GMT

前言:

《机器学习实战》笔记 (第八章和第九章) (附Python3版代码)(Machine Learning in Action)

《机器学习实战》(Machine Learning in Action)是一本常见的机器学习入门书,书中代码由Python2写成。由于现时Python2已逐渐退出舞台,所以这篇文章将该书的所有代码部分用Python3重写。
代码上传GitHub: https://github.com/kungbob/Machine_Learning_In_Action
原版Python2代码:https://www.manning.com/books/machine-learning-in-action

注: 第八章和第九章主要内容是用不同的回归方法来预测数据,可视为第五章逻辑回归的扩展。

第八章: 预测数值型数据:回归

  1. 线性回归
    优点: 结果已于理解,计算上不复杂。
    缺点: 对非线性的数据拟合不好。
    适用数据类型: 数值型和标称型数据。

一般方法

  1. 收集数据: 采用任意方法收集数据。
  2. 准备数据: 回归需要数值型数据,标称型数据将被转成二值型数据(0或1)。
  3. 分析数据: 会出数据的可视化二维图讲有助于对数值做出理解和分析,采用缩减发求得新回归系数后,可以讲拟合线在图上左对比。
  4. 训练算法: 找到回归系数。
  5. 测试算法: 使用\(R^2\)或者预测值和数据的拟合度,分析模型的效果。
  6. 使用算法: 使用回归可以在给定输入的时候预测出一个数值。

示例:用回归法预测乐高套装的价格
注: 书中的Google API并不可用。

  1. 收集数据: 用Google Shopping的API收集数据。
  2. 准备数据: 从返回的JSON数据中抽取价格。
  3. 分析数据: 可视化并观察数据。
  4. 训练算法: 构建不同的模型,采用逐步线性回归和直接的线性回归模型。
  5. 测试算法: 使用交叉验证来测试不同的模型,分析哪个效果最好。
  6. 使用算法: 生成数据模型。

regression.py完整代码(章节中Python命令行代码在GitHub):

from numpy import *
from time import sleep
import json, urllib.request

# Stand regression function and data import function, program 8_1
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

def standRegres(xArr, yArr):
    xMat = mat(xArr)
    yMat = mat(yArr).T
    xTx = xMat.T * xMat
    if linalg.det(xTx) == 0.0:
        print("This matrix is singular, cannot do inverse")
        return
    ws = xTx.I * (xMat.T * yMat)
    return ws

# weighting regression function, program 8_2
def lwlr(testPoint, xArr, yArr, k = 1.0):
    xMat = mat(xArr)
    yMat = mat(yArr).T
    m = shape(xMat)[0]
    # Create diagonal matrix
    weights = mat(eye((m)))
    for j in range(m):
        # Decrease the weighting value in expotential
        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

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

# Deducing abalone's age, chapter 8_3
def rssError(yArr, yHatArr):
    return ((yArr - yHatArr) ** 2).sum()

# Ridge regression, program 8_3
def ridgeRegres(xMat, yMat, lam = 0.2):
    xTx = xMat.T * xMat
    denom = xTx + eye(shape(xMat)[1]) * lam
    if linalg.det(denom) == 0.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
    # Data regularization
    xMeans = mean(xMat, 0)
    xVar = var(xMat, 0)
    xMat = (xMat - xMeans) / xVar
    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

# Not shown in the book, provided in source code
# Regularize by columns
def regularize(xMat):
    inMat = xMat.copy()
    # Calc mean then subtract it off
    inMeans = mean(inMat,0)
    # Calc variance of Xi then divide by it
    inVar = var(inMat,0)
    inMat = (inMat - inMeans)/inVar
    return inMat

# Forward stagewise regression, program 8_4
def stageWise(xArr, yArr, eps = 0.01, numIt = 100):
    xMat = mat(xArr)
    yMat = mat(yArr).T
    yMean = mean(yMat, 0)
    yMat = yMat - yMean
    xMat = regularize(xMat)
    m, n = shape(xMat)
    returnMat = zeros((numIt, n))
    ws = zeros((n, 1))
    wsTest = ws.copy()
    wsMax = ws.copy()
    for i in range(numIt):
        print(ws.T)
        lowestError = inf
        for j in range(n):
            for sign in [-1, 1]:
                wsTest = ws.copy()
                wsTest[j] += eps * sign
                yTest = xMat * wsTest
                rssE = rssError(yMat.A, yTest.A)
                if rssE < lowestError:
                    lowestError = rssE
                    wsMax = wsTest
        ws = wsMax.copy()
        returnMat[i, :] = ws.T
    return returnMat

# Retrieving shopping message function, program 8_5
# Outdated API, should return 404 Error
def searchForSet(retX, retY, setNum, yr, numPce, origPrc):
    sleep(10)
    myAPIstr = 'get from code.google.com'
    # Replace own API with myAPIstr
    searchURL = 'https://www.googleapis.com/shopping/search/v1/public/products?key=%s&country=US&q=lego+%d&alt=json' % (myAPIstr, setNum)
    print(searchURL)
    pg = urllib.request.urlopen(searchURL)
    retDict = json.loads(pg.read())
    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['price']
                # Filter imcomplete set
                if sellingPrice > origPrc * 0.5:
                    print("%d\t%d\t%d\t%f\t%f" %\
                        (yr, numPce, newFlag, origPrc, sellingPrice))
                    retX.append([yr, numPce, newFlag, origPrc])
                    retY.append(sellingPrice)
        except:
            print("problem with item %d" % i)

def setDataCollect(retX, retY):
    searchForSet(retX, retY, 8288, 2006, 800, 49.99)
    searchForSet(retX, retY, 10030, 2002, 3096, 269.99)
    searchForSet(retX, retY, 10179, 2007, 5195, 499.99)
    searchForSet(retX, retY, 10181, 2007, 3428, 199.99)
    searchForSet(retX, retY, 10189, 2008, 5922, 299.99)
    searchForSet(retX, retY, 10196, 2009, 3263, 249.99)

# Cross validation test on ridges regression, program 8_6
def crossValidation(xArr, yArr, numVal = 10):
    m = len(yArr)
    indexList = range(m)
    errorMat = zeros((numVal, 30))
    for i in range(numVal):
        # Create training set and testing set
        trainX = []
        trainY = []
        testX = []
        testY = []
        random.shuffle(indexList)
        for j in range(m):
            # Append data into training set and testing set
            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 = ridgeTest(trainX, trainY)
        for k in range(30):
            # Standardize all testing data using training data
            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)
            errorMat[i, k] = rssError(yEst.T.A, array(testY))
    meanErrors = mean(errorMat, 0)
    minMean = float(min(meanErrors))
    bestWeights = wMat[nonzero(meanErrors == minMean)]
    xMat = mat(xArr)
    yMat = amt(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))

第九章: 树回归

优点: 可以对复杂和非线性的数据建模。
缺点: 结果不易理解。
适用数据类型: 数值型和标称型数据。

树回归的一般方法

  1. 收集数据: 采用任意方法收集数据。
  2. 准备数据: 需要数值型的数据,标称型数据应该映射成二值型数据。
  3. 分析数据: 绘出数据的二维可视化显示结果,以字典方式生成树。
  4. 训练算法: 大部分时间都花费在叶节点树模型的构建上。
  5. 测试算法: 使用测试数据上的\(R^2\)值来分析模型的结果。
  6. 使用算法: 使用训练出的树做预测,预测结果还可以用来做很多事情。

示例: 利用GUI对回归树调优

  1. 收集数据: 提供的文本文件。
  2. 准备数据: 用Python解析上述文件,得到数值型数据。
  3. 分析数据: 用Tkinter构建一个GUI来展示模型和数据。
  4. 训练算法: 训练一颗回归树和一颗模型树,并于数据集一起展示出来。
  5. 测试算法: 不需要。
  6. 使用算法: 透过GUI测试不同参数的影响,帮助选择模型的类型。

regTrees.py完整代码(章节中Python命令行代码在GitHub):

from numpy import *

class treeNode():
    def __init__(self, feat, val, right, left):
        featureToSplitOn = feat
        valueOfSplit = val
        rightBranch = right
        leftBranch = left

# Implementation of CART algorithm, program 9_1
def loadDataSet(fileName):
    dataMat = []
    fr = open(fileName)
    for line in fr.readlines():
        curLine = line.strip().split('\t')
        # Change every line into float number
        fltLine = list(map(float, curLine))
        dataMat.append(fltLine)
    return dataMat

def binSplitDataSet(dataSet, feature, value):
    mat0 = dataSet[nonzero(dataSet[:, feature] > value)[0], :]
    mat1 = dataSet[nonzero(dataSet[:, feature] <= value)[0], :]
    return mat0, mat1

# Split function of regression Tree, program 9_2
def regLeaf(dataSet):
    return mean(dataSet[:, -1])

def regErr(dataSet):
    return var(dataSet[:, -1]) * shape(dataSet)[0]

def createTree(dataSet, leafType = regLeaf, errType = regErr, ops = (1, 4)):
    feat, val = chooseBestSplit(dataSet, leafType, errType, ops)
    # Return node value if requirements are satisfied
    if feat == None:
        return val
    retTree = {}
    retTree['spInd'] = feat
    retTree['spVal'] = val
    lSet, rSet = binSplitDataSet(dataSet, feat, val)
    retTree['left'] = createTree(lSet, leafType, errType, ops)
    retTree['right'] = createTree(rSet, leafType, errType, ops)
    return retTree

def chooseBestSplit(dataSet, leafType = regLeaf, errType = regErr, ops= (1, 4)):
    tolS = ops[0]
    tolN = ops[1]
    # If all values are equal, return
    if len(set(dataSet[:, -1].T.tolist()[0])) == 1:
        return None, leafType(dataSet)
    m, n = shape(dataSet)
    S = errType(dataSet)
    bestS = inf
    bestIndex = 0
    bestValue = 0
    for featIndex in range(n - 1):
        for splitVal in set(dataSet[:, featIndex].T.tolist()[0]):
            mat0, mat1 = binSplitDataSet(dataSet, featIndex, splitVal)
            if (shape(mat0)[0] < tolN) or (shape(mat1)[0] < tolN):
                continue
            newS = errType(mat0) + errType(mat1)
            if newS < bestS:
                bestIndex = featIndex
                bestValue = splitVal
                bestS = newS
    # If error does not reduce, return
    if (S - bestS) < tolS:
        return None, leafType(dataSet)
    mat0, mat1 = binSplitDataSet(dataSet, bestIndex, bestValue)
    # If new data set is small, return
    if (shape(mat0)[0] < tolN) or (shape(mat1)[0] < tolN):
        return None, leafType(dataSet)
    return bestIndex, bestValue

def isTree(obj):
    return (type(obj).__name__ == 'dict')

def getMean(tree):
    if isTree(tree['right']):
        tree['right'] = getMean(tree['right'])
    if isTree(tree['left']):
        tree['left'] = getMean(tree['left'])
    return (tree['left'] + tree['right']) / 2.0

def prune(tree, testData):
    # If no testing data, delete tree
    if shape(testData)[0] == 0:
        return getMean(tree)
    if (isTree(tree['right']) or isTree(tree['left'])):
        lSet, rSet = binSplitDataSet(testData, tree['spInd'], tree['spVal'])
    if isTree(tree['left']):
        tree['left'] = prune(tree['left'], lSet)
    if isTree(tree['right']):
        tree['right'] = prune(tree['right'], rSet)
    if not isTree(tree['left']) and not isTree(tree['right']):
        lSet, rSet = binSplitDataSet(testData, tree['spInd'], tree['spVal'])
        errorNoMerge = sum(power(lSet[:, -1] - tree['left'], 2)) + \
            sum(power(rSet[:, -1] - tree['right'], 2))
        treeMean = (tree['left'] + tree['right']) / 2.0
        errorMerge = sum(power(testData[:, -1] - treeMean, 2))
        if errorMerge < errorNoMerge:
            print("merging")
            return treeMean
        else:
            return treeMean
    else:
        return tree

# Model Tree's leaf node generation, program 9_4
def linearSolve(dataSet):
    m, n = shape(dataSet)
    # Set the data of X and Y
    X = mat(ones((m, n)))
    Y = mat(ones((m, 1)))
    X[:, 1:n] = dataSet[:, 0: n-1]
    Y = dataSet[:, -1]
    xTx = X.T * X
    if linalg.det(xTx) == 0.0:
        raise NameError('This matrix is singular, cannot do inverse, \n\
        try increasing the second value of ops')
    ws = xTx.I * (X.T * Y)
    return ws, X, Y

def modelLeaf(dataSet):
    ws, X, Y = linearSolve(dataSet)
    yHat = X * ws
    return sum(power(Y - yHat, 2))

def modelErr(dataSet):
    ws, X, Y = linearSolve(dataSet)
    yHat = X * ws
    return sum(power(Y - yHat, 2))

# Use tree regression to predict, program 9_5
def regTreeEval(model, inDat):
    return float(model)

def modelTreeEval(model, inDat):
    n = inDat.shape[1]
    X = mat(ones((1, n)))
    X[:, 1: n] = inDat[:, :-1]
    return float(X * model)

def treeForeCast(tree, inData, modelEval= regTreeEval):
    if not isTree(tree):
        return modelEval(tree, inData)
    if inData[tree['spInd']] > tree['spVal']:
        if isTree(tree['left']):
            return treeForeCast(tree['left'], inData, modelEval)
        else:
            return modelEval(tree['left'], inData)
    else:
        if isTree(tree['right']):
            return treeForeCast(tree['right'], inData, modelEval)
        else:
            return modelEval(tree['right'], inData)

def createForeCast(tree, testData, modelEval= regTreeEval):
    m = len(testData)
    yHat = mat(zeros((m, 1)))
    for i in range(m):
        yHat[i, 0] = treeForeCast(tree, mat(testData[i]), modelEval)
    return yHat

treeExplore.py完整代码(章节中Python命令行代码在GitHub):

from numpy import *
from tkinter import *
import regTrees
# Combine Matplotlib and tkinter, program 9_7
import matplotlib
matplotlib.use('TkAgg')
from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg
from matplotlib.figure import Figure

def reDraw(tolS, tolN):
    reDraw.f.clf()
    reDraw.a = reDraw.f.add_subplot(111)
    # Check which is selected
    if chkBtnVar.get():
        if tolN < 2:
            tolN = 2
        myTree = regTrees.createTree(reDraw.rawDat, regTrees.modelLeaf, \
                                    regTrees.modelErr, (tolS, tolN))
        yHat = regTrees.createForeCast(myTree, reDraw.testDat, \
                                        regTrees.modelTreeEval)
    else:
        myTree = regTrees.createTree(reDraw.rawDat, ops = (tolS, tolN))
        yHat = regTrees.createForeCast(myTree, reDraw.testDat)
    reDraw.a.scatter(reDraw.rawDat[:, 0].A, reDraw.rawDat[:, 1].A, s = 5)
    reDraw.a.plot(reDraw.testDat, yHat, linewidth = 2.0)
    reDraw.canvas.show()

def getInputs():
    try:
        tolN = int(tolNentry.get())
    except:
        tolN = 10
        print("enter Integer for tolN")
        # Clear wrong input and replace with default value
        tolNentry.delete(0, END)
        tolNentry.insert(0, '10')
    try:
        tolS = float(tolSentry.get())
    except:
        tolS = 1.0
        print("enter Float for tolS")
        tolSentry.delete(0, end)
        TOLsENTRY.INSERT(0, '1.0')
    return tolN, tolS

def drawNewTree():
    tolN, tolS = getInputs()
    reDraw(tolS, tolN)

root = Tk()

# Abandon in chapter 9_7_2
# Label(root, text= "Plot Place Holder").grid(row= 0, columnspan= 3)

# Added in chapter 9_7_2
reDraw.f = Figure(figsize = (5, 1), dpi = 100)
reDraw.canvas = FigureCanvasTkAgg(reDraw.f, master = root)
reDraw.canvas.show()
reDraw.canvas.get_tk_widget().grid(row = 0, columnspan = 3)

Label(root, text= "tolN").grid(row = 1, column = 0)
tolNentry = Entry(root)
tolNentry.grid(row = 1, column = 1)
tolNentry.insert(0, '10')
Label(root, text = "tolS").grid(row = 2, column = 0)
tolSentry = Entry(root)
tolSentry.grid(row = 2, column = 1)
tolSentry.insert(0, '1.0')
Button(root, text = 'ReDraw', command = drawNewTree).grid(row = 1, column = 2, \
                                                            rowspan = 3)
chkBtnVar = IntVar()
chkBtn = Checkbutton(root, text = "Model Tree", variable = chkBtnVar)
chkBtn.grid(row = 3, column = 0, columnspan = 2)

reDraw.rawDat = mat(regTrees.loadDataSet('sine.txt'))
reDraw.testDat = arange(min(reDraw.rawDat[:, 0]), \
                            max(reDraw.rawDat[:, 0]), 0.01)

reDraw(1.0, 10)
Button(root, text = 'Quit', fg = 'black', command = root.quit).grid(row = 1,
                                                                    column = 2)
root.mainloop()

]]>
<![CDATA[《机器学习实战》笔记 (第五章至第七章) (附Python3版代码)(Machine Learning in Action)]]>https://codeproducers.com/20180401/5e0456874b43c3452b5136ddSun, 01 Apr 2018 13:33:54 GMT

前言:

《机器学习实战》笔记 (第五章至第七章) (附Python3版代码)(Machine Learning in Action)

《机器学习实战》(Machine Learning in Action)是一本常见的机器学习入门书,书中代码由Python2写成。由于现时Python2已逐渐退出舞台,所以这篇文章将该书的所有代码部分用Python3重写。
代码上传GitHub: https://github.com/kungbob/Machine_Learning_In_Action
原版Python2代码:https://www.manning.com/books/machine-learning-in-action

第五章: Logistic回归

优点: 计算代价不高,易于理解和实现。
缺点: 容易欠拟合,分类精度可能不高。
适用数据类型: 数值型和标称型数据。

Logistic回归的一般过程

  1. 收集数据: 采用任意方法收集数据。
  2. 准备数据: 由于需要进行距离计算,因此要求数据类型为数值型。另外,结构化数据格式则最佳。
  3. 分析数据: 采用任意方法对数据进行分析。
  4. 训练算法: 大部分时间将用于训练,训练的目的是为了找到最佳的分类回归系统。
  5. 测试算法: 一旦训练步骤完成,分类将会很快。
  6. 使用算法: 首先需要输入一些数据,并将其转换成对应的结构化数据。接着,基于训练好的回归系数就可以对这些数值进行简单的回归计算,判定他们属于哪个类别;在这之后,就可以在输出的类型上做一些其他分析工作。

梯度下降算法:
梯度下降算法和梯度上升算法是一样的,只是公式中的加法需要变成减法。因此,梯度下降算法的公式为:
\[w := w - \alpha \nabla_w f(w) \]
梯度上升算法用来求函数的最大值,而梯度下降算法用来求函数的最小值。

示例: 使用Logistic回归估计马疝病的死亡率

  1. 收集数据: 给定数据文件。
  2. 准备数据: 用Python解析文本文件并填充缺失值。
  3. 分析数据: 可视化并观察数据。
  4. 训练算法: 使用优化算法,找到最佳的系数。
  5. 测试算法: 为了量化回归的效果,需要观测错误率。根据错误率决定是否回退到训练阶段,通过改变递送的次数和步长等参数来得到更好的回归函数。
  6. 使用算法: 实现一个简单的命令行程序来手机马的症状并输出预测结果并非难事,可以作为留给读者的习题。

logRegres.py完整代码(章节中Python命令行代码在GitHub):

from numpy import *
import matplotlib.pyplot as plt

# Logistic Regression Gradient Ascent Method, program 5_1
def loadDataSet():
    dataMat = []
    labelMat = []
    fr = open('testSet.txt')
    for line in fr.readlines():
        lineArr = line.strip().split()
        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
        labelMat.append(int(lineArr[2]))
    return dataMat, labelMat

def sigmoid(inX):
    return 1.0 / (1 + exp(-inX))

def gradAscent(dataMatIn, classLabels):
    # Convert the matrix into Numpy format
    dataMatrix = mat(dataMatIn)
    labelMat = mat(classLabels).transpose()
    m, n = shape(dataMatrix)
    alpha = 0.001
    maxCycles = 500
    weights = ones((n, 1))
    # Matrix multiplication
    for k in range(maxCycles):
        h = sigmoid(dataMatrix * weights)
        error = (labelMat - h)
        weights = weights + alpha * dataMatrix.transpose() * error
    return weights

# Ploting the dataset and the decision boundary
def plotBestFit(weights):
    dataMat, labelMat = loadDataSet()
    dataArr = array(dataMat)
    n = shape(dataArr)[0]
    xcord1 = []
    ycord1 = []
    xcord2 = []
    ycord2 = []
    for i in range(n):
        if int(labelMat[i]) == i:
            xcord1.append(dataArr[i, 1])
            ycord1.append(dataArr[i, 2])
        else:
            xcord2.append(dataArr[i, 1])
            ycord2.append(dataArr[i, 2])
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(xcord1, ycord1, s = 30, c = 'red', marker = 's')
    ax.scatter(xcord2, ycord2, s = 30, c = 'green')
    x = arange(-3.0, 3.0, 0.1)
    # Best fiting line
    y = (-weights[0] - weights[1] * x) / weights[2]
    ax.plot(x, y)
    plt.xlabel('X1')
    plt.ylabel('X2')
    plt.show()

# Stochastic gradient ascent, program 5_3
def stocGradAscent0(dataMatrix, classLabels):
    m, n = shape(dataMatrix)
    alpha = 0.01
    weights = ones(n)
    for i in range(m):
        h = sigmoid(sum(dataMatrix[i] * weights))
        error = classLabels[i] - h
        weights = weights + alpha * error * dataMatrix[i]
    return weights

# Improved Stochastic gradient ascent, program 5_4
def stocGradAscent1(dataMatrix, classLabels, numIter = 150):
    m, n = shape(dataMatrix)
    weights = ones(n)
    for j in range(numIter):
        dataIndex = list(range(m))
        for i in range(m):
            # Modifying alpha for each iteration
            alpha = 4 / (1.0 + j + i) + 0.01
            # Randomly pick a data to update
            randIndex = int(random.uniform(0, len(dataIndex)))
            h = sigmoid(sum(dataMatrix[randIndex] * weights))
            error = classLabels[randIndex] - h
            weights = weights + alpha * error * dataMatrix[randIndex]
            del(dataIndex[randIndex])
    return weights

# Logistic regession classifier, program 5_5
def classifyVector(inX, weights):
    prob = sigmoid(sum(inX * weights))
    if prob > 0.5:
        return 1.0
    else:
        return 0.0

def colicTest():
    frTrain = open('horseColicTraining.txt')
    frTest = open('horseColicTest.txt')
    trainingSet = []
    trainingLabels = []
    for line in frTrain.readlines():
        currLine = line.strip().split('\t')
        lineArr = []
        for i in range(21):
            lineArr.append(float(currLine[i]))
        trainingSet.append(lineArr)
        trainingLabels.append(float(currLine[21]))
    trainWeights = stocGradAscent1(array(trainingSet), trainingLabels, 500)
    errorCount = 0
    numTestVec = 0.0
    for line in frTest.readlines():
        numTestVec += 1.0
        currLine = line.strip().split('\t')
        lineArr = []
        for i in range(21):
            lineArr.append(float(currLine[i]))
        if int(classifyVector(array(lineArr), trainWeights)) != int(currLine[21]):
            errorCount += 1
    errorRate = (float(errorCount) / numTestVec)
    print("the error rate of this test is: %f" % errorRate)
    return errorRate

def multiTest():
    numTests = 10
    errorSum = 0.0
    for k in range(numTests):
        errorSum += colicTest()
    print("after %d iterations the average error rate is: %f" % (numTests, errorSum/float(numTests)))

第六章: 支持向量机

优点: 泛化错误低(Generalization error),计算开销不大,结果易解释。
缺点: 对参数调节和核函数的选择敏感,原始分类器不加修改仅使用于处理二类问题(Two-class Classification)。
适用数据类型: 数值型和标称型数据。

SVM的一般流程:

  1. 收集数据: 可以使用任意方法。
  2. 准备数据: 需要数值型数据。
  3. 分析数据: 有助于可视化分隔超平面。
  4. 训练算法: SVM的大部分时间都源自训练,该过程主要实现两个参数的调优。
  5. 测试算法: 十分简单的计算过程就可以实现。
  6. 使用算法: 几乎所有分类问题都可以使用SVM,值得一提的是,SVM本身是一个二类分类器,对多类问题应用SVM需要对代码做一些修改。

示例:基于SVM的数字识别

  1. 收集数据: 提供的文本文件。
  2. 准备数据: 基于二值图像构造向量。
  3. 分析数据: 对图像向量进行目测。
  4. 训练算法: 采用两种不同的核函数,并对径向核函数采用不同的设置来运行SMO算法。
  5. 测试算法: 编写一个函数来测试不同的核函数并计算错误率。
  6. 使用算法: 一个图像识别的完整应用还需要一些图像处理的知识,这里并不深入介绍。

svmMLiA.py完整代码(章节中Python命令行代码在GitHub):

from numpy import *
from os import listdir

# Supporting functions in SMO algorithm, program 6_1
def loadDataSet(fileName):
    dataMat = []
    labelMat = []
    fr = open(fileName)
    for line in fr.readlines():
        lineArr = line.strip().split('\t')
        dataMat.append([float(lineArr[0]), float(lineArr[1])])
        labelMat.append(float(lineArr[2]))
    return dataMat, labelMat

def selectJrand(i, m):
    j = i
    while (j == i):
        j = int(random.uniform(0, m))
    return j

def clipAlpha(aj, H, L):
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj

# simplified SMO algorithm, program 6_2
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
    dataMatrix = mat(dataMatIn)
    labelMat = mat(classLabels).transpose()
    b = 0
    m, n = shape(dataMatrix)
    alphas = mat(zeros((m, 1)))
    iter = 0
    while (iter < maxIter):
        alphaPairsChanged = 0
        for i in range(m):
            fXi = float(multiply(alphas, labelMat).T*\
                        (dataMatrix*dataMatrix[i, :].T)) + b
            Ei = fXi - float(labelMat[i])
            # Test the value of alpha whether it needs optimization.
            if ((labelMat[i] * Ei < -toler) and (alphas[i] < C)) or \
                ((labelMat[i] * Ei > toler) and \
                (alphas[i] > 0)):
                # Randomly pick the second alpha
                j = selectJrand(i, m)
                fXj = float(multiply(alphas, labelMat).T*\
                            (dataMatrix*dataMatrix[j, :].T)) + b
                Ej = fXj - float(labelMat[j])
                alphaIold = alphas[i].copy()
                alphaJold = alphas[j].copy()
                # Ensure the range of alpha is between 0 to C
                if (labelMat[i] != labelMat[j]):
                    L = max(0, alphas[j] - alphas[i])
                    H = min(C, C + alphas[j] - alphas[i])
                else:
                    L = max(0, alphas[j] + alphas[i] - C)
                    H = min(C, alphas[j] + alphas[i])
                if L == H:
                    print("L == H")
                    continue
                eta = 2.0 * dataMatrix[i, :] * dataMatrix[j, :].T - \
                      dataMatrix[i, :] * dataMatrix[i, :].T - \
                      dataMatrix[j, :] * dataMatrix[j, :].T
                if eta >= 0:
                    print("eta>=0")
                    continue
                alphas[j] -= labelMat[j] * (Ei - Ej) / eta
                alphas[j] = clipAlpha(alphas[j], H, L)
                if (abs(alphas[j] - alphaJold) < 0.0001):
                    print("j not moving enough")
                    continue
                # Modify i with value equals to j, but in opposite direction
                alphas[i] += labelMat[j] * labelMat[i] * \
                            (alphaJold - alphas[j])
                b1 = b - Ei - labelMat[i] * (alphas[i] - alphaIold) * \
                    dataMatrix[i, :] * dataMatrix[i, :].T - \
                    labelMat[j] * (alphas[j] - alphaJold) * \
                    dataMatrix[i, :] * dataMatrix[j, :].T
                # Set constant term
                b2 = b - Ej - labelMat[i] * (alphas[i] - alphaIold) * \
                     dataMatrix[i, :] * dataMatrix[j, :].T - \
                     labelMat[j] * (alphas[j] - alphaJold) * \
                     dataMatrix[j, :] * dataMatrix[j, :].T
                if (0 < alphas[i]) and (C > alphas[i]):
                    b = b1
                elif (0 < alphas[j]) and (C > alphas[j]):
                    b = b2
                else:
                    b = (b1 + b2) / 2.0
                alphaPairsChanged += 1
                print("iter: %d i: %d, pairs changed %d" % \
                        (iter, i, alphaPairsChanged))
        if (alphaPairsChanged == 0):
            iter += 1
        else:
            iter = 0
        print("iteration number: %d" % iter)
    return b, alphas

# Complete version of Platt SMO Supporting functions, program 6_3
'''
class optStruct:
    def __init__(self, dataMatIn, classLabels, C, toler):
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = shape(dataMatIn)[0]
        self.alphas = mat(zeros((self.m, 1)))
        self.b = 0
        # Save the error into Cache
        self.eCache = mat(zeros((self.m, 2)))
'''

# Modified for Kernal function, program 6_6
class optStruct:
    def __init__(self, dataMatIn, classLabels, C, toler, kTup):
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = shape(dataMatIn)[0]
        self.alphas = mat(zeros((self.m, 1)))
        self.b = 0
        # Save the error into Cache
        self.eCache = mat(zeros((self.m, 2)))
        self.K = mat(zeros((self.m, self.m)))
        for i in range(self.m):
            self.K[:, i] = n (self.X, self.X[i, :], kTup)

'''
def calcEk(oS, k):
    fXk = float(multiply(oS.alphas, oS.labelMat).T*\
            (oS.X*oS.X[k, :].T)) + oS.b
    Ek = fXk - float(oS.labelMat[k])
    return Ek
'''

# Modified for kernal function, program 6_7
def calcEk(oS, k):
    fXk = float(multiply(oS.alphas, oS.labelMat).T * oS.K[:, k] + oS.b)
    Ek = fXk - float(oS.labelMat[k])
    return Ek

# Preparation for the inner interation loop
def selectJ(i, oS, Ei):
    maxK = -1
    maxDeltaE = 0
    Ej = 0
    oS.eCache[i] = [1, Ei]
    validEcacheList = nonzero(oS.eCache[:, 0].A)[0]
    if (len(validEcacheList)) > 1:
        for k in validEcacheList:
            if k == i:
                continue
            Ek = calcEk(oS, k)
            deltaE = abs(Ei - Ek)
            # Select the largest step distance of j
            if (deltaE > maxDeltaE):
                maxK = k
                maxDeltaE = deltaE
                Ej = Ek
        return maxK, Ej
    else:
        j = selectJrand(i, oS.m)
        Ej = calcEk(oS, j)
    return j, Ej

def updateEk(oS, k):
    Ek = calcEk(oS, k)
    oS.eCache[k] = [1, Ek]

# Complete Platt SMO algorithm, optimised version, program 6_4
'''
def innerL(i, oS):
    Ei = calcEk(oS, i)
    if ((oS.labelMat[i]*Ei < -oS.tol) and (oS.alphas[i] < oS.C)) or \
        ((oS.labelMat[i]*Ei > oS.tol) and (oS.alphas[i] > 0)):
        j, Ej = selectJ(i, oS, Ei)
        # Select the second alpha
        alphaIold = oS.alphas[i].copy()
        alphaJold = oS.alphas[j].copy()
        if (oS.labelMat[i] != oS.labelMat[j]):
            L = max(0, oS.alphas[j] - oS.alphas[i])
            H = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i])
        else:
            L = max(0, oS.alphas[j] + oS.alphas[i] - oS.C)
            H = min(oS.C, oS.alphas[j] + oS.alphas[i])
        if L == H:
            print("L == H")
            return 0
        eta = 2.0 * oS.X[i, :] * oS.X[j, :].T -oS.X[i, :] * oS.X[i, :].T - \
                oS.X[j, :] * oS.X[j, :].T
        if eta >= 0:
            print("eta >= 0")
            return 0
        oS.alphas[j] -= oS.labelMat[j]* (Ei - Ej) / eta
        oS.alphas[j] = clipAlpha(oS.alphas[j], H, L)
        updateEk(oS, j)
        if (abs(oS.alphas[j] - alphaJold) < 0.0001):
            print("j not moving enough")
            return 0
        oS.alphas[i] += oS.labelMat[j] * oS.labelMat[i] * \
                        (alphaJold - oS.alphas[j])
        updateEk(oS, i)
        b1 = oS.b - Ei - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * \
                oS.X[i, :] * oS.X[i, :].T - oS.labelMat[j] * \
                (oS.alphas[j] - alphaJold) * oS.X[i, :] * oS.X[j, :].T
        b2 = oS.b - Ej - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * \
                oS.X[i, :] * oS.X[j, :].T - oS.labelMat[j] * \
                (oS.alphas[j] - alphaJold) * oS.X[j, :] * oS.X[j, :].T
        if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]):
            oS.b = b1
        elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]) :
            oS.b = b2
        else:
            oS.b = (b1 + b2) / 2.0
        return 1
    else:
        return 0
'''

# Inner iteration for Kernal function, program 6_7
def innerL(i, oS):
    Ei = calcEk(oS, i)
    if ((oS.labelMat[i]*Ei < -oS.tol) and (oS.alphas[i] < oS.C)) or \
        ((oS.labelMat[i]*Ei > oS.tol) and (oS.alphas[i] > 0)):
        j, Ej = selectJ(i, oS, Ei)
        # Select the second alpha
        alphaIold = oS.alphas[i].copy()
        alphaJold = oS.alphas[j].copy()
        if (oS.labelMat[i] != oS.labelMat[j]):
            L = max(0, oS.alphas[j] - oS.alphas[i])
            H = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i])
        else:
            L = max(0, oS.alphas[j] + oS.alphas[i] - oS.C)
            H = min(oS.C, oS.alphas[j] + oS.alphas[i])
        if L == H:
            print("L == H")
            return 0
        # Modified in program 6_7
        eta = 2.0 * oS.K[i, j] - oS.K[i, i] - oS.K[j, j]
        if eta >= 0:
            print("eta >= 0")
            return 0
        oS.alphas[j] -= oS.labelMat[j]* (Ei - Ej) / eta
        oS.alphas[j] = clipAlpha(oS.alphas[j], H, L)
        updateEk(oS, j)
        if (abs(oS.alphas[j] - alphaJold) < 0.0001):
            print("j not moving enough")
            return 0
        oS.alphas[i] += oS.labelMat[j] * oS.labelMat[i] * \
                        (alphaJold - oS.alphas[j])
        updateEk(oS, i)
        # Modified in program 6_7
        b1 = oS.b - Ei - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * \
                oS.K[i, i] - oS.labelMat[j] * (oS.alphas[j] - alphaJold) * \
                oS.K[i, j]
        b2 = oS.b - Ej - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * \
                oS.K[i, j] - oS.labelMat[j] * (oS.alphas[j] - alphaJold) * \
                oS.K[j, j]

        if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]):
            oS.b = b1
        elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]) :
            oS.b = b2
        else:
            oS.b = (b1 + b2) / 2.0
        return 1
    else:
        return 0

# Complete Platt SMO outer iteration, program 6_5
def smoP(dataMatIn, classLabels, C, toler, maxIter, kTup=('lin', 0)):
    # Add kTup parameter at the end
    oS = optStruct(mat(dataMatIn), mat(classLabels).transpose(), C, toler, kTup)
    iter = 0
    entireSet = True
    alphaPairsChanged = 0
    # Iterate all value
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
        alphaPairsChanged = 0
        if entireSet:
            for i in range(oS.m):
                alphaPairsChanged += innerL(i, oS)
                print("fullSet, iter: %d i:%d, pairs changed %d" %\
                        (iter, i, alphaPairsChanged))
            iter += 1
        else:
            # Iterate all non-boundaries value
            nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]
            for i in nonBoundIs:
                alphaPairsChanged += innerL(i, oS)
                print("non-bound, iter: %d i: %d, pairs changed %d" % \
                        (iter, i, alphaPairsChanged))
            iter += 1
        if entireSet:
            entireSet = False
        elif (alphaPairsChanged == 0):
            entireSet = True
        print("iteration number: %d" % iter)
    return oS.b, oS.alphas

# Calculate w value, chapter 6_4
def calcWs(alphas, dataArr, classLabels):
    X = mat(dataArr)
    labelMat = mat(classLabels).transpose()
    m,n = shape(X)
    w = zeros((n, 1))
    for i in range(m):
        w += multiply(alphas[i] * labelMat[i], X[i, :].T)
    return w

# Kernal translation function, program 6_6
def kernelTrans(X, A, kTup):
    m, n = shape(X)
    K = mat(zeros((m, 1)))
    if kTup[0] == 'lin':
        K = X * A.T
    elif kTup[0] == 'rbf':
        for j in range(m):
            deltaRow = X[j, :] - A
            K[j] = deltaRow * deltaRow.T
        # Division betwen elements
        K = exp(K /(-1*kTup[1]**2))
    else:
        raise NameError('Houston We Have a Problem -- \
        That Kernal is not recongnized')
    return K

def testRbf(k1 = 1.3):
    dataArr, labelArr = loadDataSet('testSetRBF.txt')
    b, alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', k1))
    datMat = mat(dataArr)
    labelMat = mat(labelArr).transpose()
    svInd = nonzero(alphas.A > 0)[0]
    # Construct SVM matrix
    sVs = datMat[svInd]
    labelSV = labelMat[svInd]
    print("there are %d Support Vectors" % shape(sVs)[0])
    m, n = shape(datMat)
    errorCount = 0
    for i in range(m):
        kernelEval = kernelTrans(sVs, datMat[i, :], ('rbf', k1))
        predict = kernelEval.T * multiply(labelSV, alphas[svInd]) + b
        if sign(predict) != sign(labelArr[i]):
            errorCount += 1
    print("the training error rate is: %f" % (float(errorCount)/m))
    dataArr, labelArr = loadDataSet('testSetRBF2.txt')
    errorCount = 0
    datMat = mat(dataArr)
    labelMat = mat(labelArr).transpose()
    m, n = shape(datMat)
    for i in range(m):
        kernelEval = kernelTrans(sVs, datMat[i, :], ('rbf', k1))
        predict = kernelEval.T * multiply(labelSV, alphas[svInd]) + b
        if sign(predict) != sign(labelArr[i]):
            errorCount += 1
    print("the test error rate is: %f" %(float(errorCount) / m))

# Handwriting recongnition based on SVM, program 6_9
def loadImages(dirName):
    hwLabels = []
    trainingFileList = listdir(dirName)
    m = len(trainingFileList)
    trainingMat = zeros((m, 1024))
    for i in range(m):
        fileNameStr = trainingFileList[i]
        fileStr = fileNameStr.split('.')[0]
        classNumStr = int(fileStr.split('_')[0])
        if classNumStr == 9:
            hwLabels.append(-1)
        else:
            hwLabels.append(1)
        trainingMat[i, :] = img2vector('%s/%s' % (dirName, fileNameStr))
    return trainingMat, hwLabels

def testDigits(kTup = ('rbf', 10)):
    dataArr, labelArr = loadImages('trainingDigits')
    b, alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, kTup)
    datMat = mat(dataArr)
    labelMat = mat(labelArr).transpose()
    svInd = nonzero(alphas.A > 0)[0]
    sVs = datMat[svInd]
    labelSV = labelMat[svInd]
    print("there are %d Support Vectors" % shape(sVs)[0])
    m, n = shape(datMat)
    errorCount = 0
    for i in range(m):
        kernelEval = kernelTrans(sVs, datMat[i, :], kTup)
        predict = kernelEval.T * multiply(labelSV, alphas[svInd]) + b
        if sign(predict) != sign(labelArr[i]):
            errorCount += 1
    print("the training error rate is: %f" % (float(errorCount) / m))
    dataArr, labelArr = loadImages('testDigits')
    errorCount = 0
    datMat = mat(dataArr)
    labelMat = mat(labelArr).transpose()
    m, n = shape(datMat)
    for i in range(m):
        kernelEval = kernelTrans(sVs, datMat[i, :], kTup)
        predict = kernelEval.T * multiply(labelSV, alphas[svInd]) + b
        if sign(predict) != sign(labelArr[i]):
            errorCount += 1
    print("test test error rate is : %f" % (float(errorCount) / m))

def img2vector(filename):
    returnVect = zeros((1, 1024))
    fr = open(filename)
    for i in range(32):
        lineStr = fr.readline()
        for j in range(32):
            returnVect[0, 32*i + j] = int(lineStr[j])
    return returnVect

第六章: 支持向量机

优点: 泛化错误低,计算开销不大,结果容易解释。
缺点: 对参数调节和和函数的选择敏感,原始分类器不加修改仅适用于处理二类问题。
使用数据类型: 数值型和标称型数据。

SVM的一般流程

  1. 收集数据: 可以使用任意方法。
  2. 准备数据: 需要数值型数据。
  3. 分析数据: 有助于可视化分隔超平面。
  4. 训练算法: SVM的大部分时间都源自训练,该过程主要实现两个参数的调优。
  5. 测试算法: 十分简单的计算过程就可以实现。
  6. 使用算法: 几乎所有分类问题都可以使用SVM。(SVM本身是一个二类分类器,对多类问题应用SVM需要对代码做一些修改)

注: 书中作者使用SMO算法作为SVM的训练算法。

svmMLiA.py完整代码(章节中Python命令行代码在GitHub):

from numpy import *
from os import listdir

# Supporting functions in SMO algorithm, program 6_1
def loadDataSet(fileName):
    dataMat = []
    labelMat = []
    fr = open(fileName)
    for line in fr.readlines():
        lineArr = line.strip().split('\t')
        dataMat.append([float(lineArr[0]), float(lineArr[1])])
        labelMat.append(float(lineArr[2]))
    return dataMat, labelMat

def selectJrand(i, m):
    j = i
    while (j == i):
        j = int(random.uniform(0, m))
    return j

def clipAlpha(aj, H, L):
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj

# simplified SMO algorithm, program 6_2
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
    dataMatrix = mat(dataMatIn)
    labelMat = mat(classLabels).transpose()
    b = 0
    m, n = shape(dataMatrix)
    alphas = mat(zeros((m, 1)))
    iter = 0
    while (iter < maxIter):
        alphaPairsChanged = 0
        for i in range(m):
            fXi = float(multiply(alphas, labelMat).T*\
                        (dataMatrix*dataMatrix[i, :].T)) + b
            Ei = fXi - float(labelMat[i])
            # Test the value of alpha whether it needs optimization.
            if ((labelMat[i] * Ei < -toler) and (alphas[i] < C)) or \
                ((labelMat[i] * Ei > toler) and \
                (alphas[i] > 0)):
                # Randomly pick the second alpha
                j = selectJrand(i, m)
                fXj = float(multiply(alphas, labelMat).T*\
                            (dataMatrix*dataMatrix[j, :].T)) + b
                Ej = fXj - float(labelMat[j])
                alphaIold = alphas[i].copy()
                alphaJold = alphas[j].copy()
                # Ensure the range of alpha is between 0 to C
                if (labelMat[i] != labelMat[j]):
                    L = max(0, alphas[j] - alphas[i])
                    H = min(C, C + alphas[j] - alphas[i])
                else:
                    L = max(0, alphas[j] + alphas[i] - C)
                    H = min(C, alphas[j] + alphas[i])
                if L == H:
                    print("L == H")
                    continue
                eta = 2.0 * dataMatrix[i, :] * dataMatrix[j, :].T - \
                      dataMatrix[i, :] * dataMatrix[i, :].T - \
                      dataMatrix[j, :] * dataMatrix[j, :].T
                if eta >= 0:
                    print("eta>=0")
                    continue
                alphas[j] -= labelMat[j] * (Ei - Ej) / eta
                alphas[j] = clipAlpha(alphas[j], H, L)
                if (abs(alphas[j] - alphaJold) < 0.0001):
                    print("j not moving enough")
                    continue
                # Modify i with value equals to j, but in opposite direction
                alphas[i] += labelMat[j] * labelMat[i] * \
                            (alphaJold - alphas[j])
                b1 = b - Ei - labelMat[i] * (alphas[i] - alphaIold) * \
                    dataMatrix[i, :] * dataMatrix[i, :].T - \
                    labelMat[j] * (alphas[j] - alphaJold) * \
                    dataMatrix[i, :] * dataMatrix[j, :].T
                # Set constant term
                b2 = b - Ej - labelMat[i] * (alphas[i] - alphaIold) * \
                     dataMatrix[i, :] * dataMatrix[j, :].T - \
                     labelMat[j] * (alphas[j] - alphaJold) * \
                     dataMatrix[j, :] * dataMatrix[j, :].T
                if (0 < alphas[i]) and (C > alphas[i]):
                    b = b1
                elif (0 < alphas[j]) and (C > alphas[j]):
                    b = b2
                else:
                    b = (b1 + b2) / 2.0
                alphaPairsChanged += 1
                print("iter: %d i: %d, pairs changed %d" % \
                        (iter, i, alphaPairsChanged))
        if (alphaPairsChanged == 0):
            iter += 1
        else:
            iter = 0
        print("iteration number: %d" % iter)
    return b, alphas

# Complete version of Platt SMO Supporting functions, program 6_3
'''
class optStruct:
    def __init__(self, dataMatIn, classLabels, C, toler):
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = shape(dataMatIn)[0]
        self.alphas = mat(zeros((self.m, 1)))
        self.b = 0
        # Save the error into Cache
        self.eCache = mat(zeros((self.m, 2)))
'''

# Modified for Kernal function, program 6_6
class optStruct:
    def __init__(self, dataMatIn, classLabels, C, toler, kTup):
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = shape(dataMatIn)[0]
        self.alphas = mat(zeros((self.m, 1)))
        self.b = 0
        # Save the error into Cache
        self.eCache = mat(zeros((self.m, 2)))
        self.K = mat(zeros((self.m, self.m)))
        for i in range(self.m):
            self.K[:, i] = n (self.X, self.X[i, :], kTup)

'''
def calcEk(oS, k):
    fXk = float(multiply(oS.alphas, oS.labelMat).T*\
            (oS.X*oS.X[k, :].T)) + oS.b
    Ek = fXk - float(oS.labelMat[k])
    return Ek
'''

# Modified for kernal function, program 6_7
def calcEk(oS, k):
    fXk = float(multiply(oS.alphas, oS.labelMat).T * oS.K[:, k] + oS.b)
    Ek = fXk - float(oS.labelMat[k])
    return Ek

# Preparation for the inner interation loop
def selectJ(i, oS, Ei):
    maxK = -1
    maxDeltaE = 0
    Ej = 0
    oS.eCache[i] = [1, Ei]
    validEcacheList = nonzero(oS.eCache[:, 0].A)[0]
    if (len(validEcacheList)) > 1:
        for k in validEcacheList:
            if k == i:
                continue
            Ek = calcEk(oS, k)
            deltaE = abs(Ei - Ek)
            # Select the largest step distance of j
            if (deltaE > maxDeltaE):
                maxK = k
                maxDeltaE = deltaE
                Ej = Ek
        return maxK, Ej
    else:
        j = selectJrand(i, oS.m)
        Ej = calcEk(oS, j)
    return j, Ej

def updateEk(oS, k):
    Ek = calcEk(oS, k)
    oS.eCache[k] = [1, Ek]

# Complete Platt SMO algorithm, optimised version, program 6_4
'''
def innerL(i, oS):
    Ei = calcEk(oS, i)
    if ((oS.labelMat[i]*Ei < -oS.tol) and (oS.alphas[i] < oS.C)) or \
        ((oS.labelMat[i]*Ei > oS.tol) and (oS.alphas[i] > 0)):
        j, Ej = selectJ(i, oS, Ei)
        # Select the second alpha
        alphaIold = oS.alphas[i].copy()
        alphaJold = oS.alphas[j].copy()
        if (oS.labelMat[i] != oS.labelMat[j]):
            L = max(0, oS.alphas[j] - oS.alphas[i])
            H = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i])
        else:
            L = max(0, oS.alphas[j] + oS.alphas[i] - oS.C)
            H = min(oS.C, oS.alphas[j] + oS.alphas[i])
        if L == H:
            print("L == H")
            return 0
        eta = 2.0 * oS.X[i, :] * oS.X[j, :].T -oS.X[i, :] * oS.X[i, :].T - \
                oS.X[j, :] * oS.X[j, :].T
        if eta >= 0:
            print("eta >= 0")
            return 0
        oS.alphas[j] -= oS.labelMat[j]* (Ei - Ej) / eta
        oS.alphas[j] = clipAlpha(oS.alphas[j], H, L)
        updateEk(oS, j)
        if (abs(oS.alphas[j] - alphaJold) < 0.0001):
            print("j not moving enough")
            return 0
        oS.alphas[i] += oS.labelMat[j] * oS.labelMat[i] * \
                        (alphaJold - oS.alphas[j])
        updateEk(oS, i)
        b1 = oS.b - Ei - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * \
                oS.X[i, :] * oS.X[i, :].T - oS.labelMat[j] * \
                (oS.alphas[j] - alphaJold) * oS.X[i, :] * oS.X[j, :].T
        b2 = oS.b - Ej - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * \
                oS.X[i, :] * oS.X[j, :].T - oS.labelMat[j] * \
                (oS.alphas[j] - alphaJold) * oS.X[j, :] * oS.X[j, :].T
        if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]):
            oS.b = b1
        elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]) :
            oS.b = b2
        else:
            oS.b = (b1 + b2) / 2.0
        return 1
    else:
        return 0
'''

# Inner iteration for Kernal function, program 6_7
def innerL(i, oS):
    Ei = calcEk(oS, i)
    if ((oS.labelMat[i]*Ei < -oS.tol) and (oS.alphas[i] < oS.C)) or \
        ((oS.labelMat[i]*Ei > oS.tol) and (oS.alphas[i] > 0)):
        j, Ej = selectJ(i, oS, Ei)
        # Select the second alpha
        alphaIold = oS.alphas[i].copy()
        alphaJold = oS.alphas[j].copy()
        if (oS.labelMat[i] != oS.labelMat[j]):
            L = max(0, oS.alphas[j] - oS.alphas[i])
            H = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i])
        else:
            L = max(0, oS.alphas[j] + oS.alphas[i] - oS.C)
            H = min(oS.C, oS.alphas[j] + oS.alphas[i])
        if L == H:
            print("L == H")
            return 0
        # Modified in program 6_7
        eta = 2.0 * oS.K[i, j] - oS.K[i, i] - oS.K[j, j]
        if eta >= 0:
            print("eta >= 0")
            return 0
        oS.alphas[j] -= oS.labelMat[j]* (Ei - Ej) / eta
        oS.alphas[j] = clipAlpha(oS.alphas[j], H, L)
        updateEk(oS, j)
        if (abs(oS.alphas[j] - alphaJold) < 0.0001):
            print("j not moving enough")
            return 0
        oS.alphas[i] += oS.labelMat[j] * oS.labelMat[i] * \
                        (alphaJold - oS.alphas[j])
        updateEk(oS, i)
        # Modified in program 6_7
        b1 = oS.b - Ei - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * \
                oS.K[i, i] - oS.labelMat[j] * (oS.alphas[j] - alphaJold) * \
                oS.K[i, j]
        b2 = oS.b - Ej - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * \
                oS.K[i, j] - oS.labelMat[j] * (oS.alphas[j] - alphaJold) * \
                oS.K[j, j]

        if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]):
            oS.b = b1
        elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]) :
            oS.b = b2
        else:
            oS.b = (b1 + b2) / 2.0
        return 1
    else:
        return 0

# Complete Platt SMO outer iteration, program 6_5
def smoP(dataMatIn, classLabels, C, toler, maxIter, kTup=('lin', 0)):
    # Add kTup parameter at the end
    oS = optStruct(mat(dataMatIn), mat(classLabels).transpose(), C, toler, kTup)
    iter = 0
    entireSet = True
    alphaPairsChanged = 0
    # Iterate all value
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
        alphaPairsChanged = 0
        if entireSet:
            for i in range(oS.m):
                alphaPairsChanged += innerL(i, oS)
                print("fullSet, iter: %d i:%d, pairs changed %d" %\
                        (iter, i, alphaPairsChanged))
            iter += 1
        else:
            # Iterate all non-boundaries value
            nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]
            for i in nonBoundIs:
                alphaPairsChanged += innerL(i, oS)
                print("non-bound, iter: %d i: %d, pairs changed %d" % \
                        (iter, i, alphaPairsChanged))
            iter += 1
        if entireSet:
            entireSet = False
        elif (alphaPairsChanged == 0):
            entireSet = True
        print("iteration number: %d" % iter)
    return oS.b, oS.alphas

# Calculate w value, chapter 6_4
def calcWs(alphas, dataArr, classLabels):
    X = mat(dataArr)
    labelMat = mat(classLabels).transpose()
    m,n = shape(X)
    w = zeros((n, 1))
    for i in range(m):
        w += multiply(alphas[i] * labelMat[i], X[i, :].T)
    return w

# Kernal translation function, program 6_6
def kernelTrans(X, A, kTup):
    m, n = shape(X)
    K = mat(zeros((m, 1)))
    if kTup[0] == 'lin':
        K = X * A.T
    elif kTup[0] == 'rbf':
        for j in range(m):
            deltaRow = X[j, :] - A
            K[j] = deltaRow * deltaRow.T
        # Division betwen elements
        K = exp(K /(-1*kTup[1]**2))
    else:
        raise NameError('Houston We Have a Problem -- \
        That Kernal is not recongnized')
    return K

def testRbf(k1 = 1.3):
    dataArr, labelArr = loadDataSet('testSetRBF.txt')
    b, alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', k1))
    datMat = mat(dataArr)
    labelMat = mat(labelArr).transpose()
    svInd = nonzero(alphas.A > 0)[0]
    # Construct SVM matrix
    sVs = datMat[svInd]
    labelSV = labelMat[svInd]
    print("there are %d Support Vectors" % shape(sVs)[0])
    m, n = shape(datMat)
    errorCount = 0
    for i in range(m):
        kernelEval = kernelTrans(sVs, datMat[i, :], ('rbf', k1))
        predict = kernelEval.T * multiply(labelSV, alphas[svInd]) + b
        if sign(predict) != sign(labelArr[i]):
            errorCount += 1
    print("the training error rate is: %f" % (float(errorCount)/m))
    dataArr, labelArr = loadDataSet('testSetRBF2.txt')
    errorCount = 0
    datMat = mat(dataArr)
    labelMat = mat(labelArr).transpose()
    m, n = shape(datMat)
    for i in range(m):
        kernelEval = kernelTrans(sVs, datMat[i, :], ('rbf', k1))
        predict = kernelEval.T * multiply(labelSV, alphas[svInd]) + b
        if sign(predict) != sign(labelArr[i]):
            errorCount += 1
    print("the test error rate is: %f" %(float(errorCount) / m))

# Handwriting recongnition based on SVM, program 6_9
def loadImages(dirName):
    hwLabels = []
    trainingFileList = listdir(dirName)
    m = len(trainingFileList)
    trainingMat = zeros((m, 1024))
    for i in range(m):
        fileNameStr = trainingFileList[i]
        fileStr = fileNameStr.split('.')[0]
        classNumStr = int(fileStr.split('_')[0])
        if classNumStr == 9:
            hwLabels.append(-1)
        else:
            hwLabels.append(1)
        trainingMat[i, :] = img2vector('%s/%s' % (dirName, fileNameStr))
    return trainingMat, hwLabels

def testDigits(kTup = ('rbf', 10)):
    dataArr, labelArr = loadImages('trainingDigits')
    b, alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, kTup)
    datMat = mat(dataArr)
    labelMat = mat(labelArr).transpose()
    svInd = nonzero(alphas.A > 0)[0]
    sVs = datMat[svInd]
    labelSV = labelMat[svInd]
    print("there are %d Support Vectors" % shape(sVs)[0])
    m, n = shape(datMat)
    errorCount = 0
    for i in range(m):
        kernelEval = kernelTrans(sVs, datMat[i, :], kTup)
        predict = kernelEval.T * multiply(labelSV, alphas[svInd]) + b
        if sign(predict) != sign(labelArr[i]):
            errorCount += 1
    print("the training error rate is: %f" % (float(errorCount) / m))
    dataArr, labelArr = loadImages('testDigits')
    errorCount = 0
    datMat = mat(dataArr)
    labelMat = mat(labelArr).transpose()
    m, n = shape(datMat)
    for i in range(m):
        kernelEval = kernelTrans(sVs, datMat[i, :], kTup)
        predict = kernelEval.T * multiply(labelSV, alphas[svInd]) + b
        if sign(predict) != sign(labelArr[i]):
            errorCount += 1
    print("test test error rate is : %f" % (float(errorCount) / m))

def img2vector(filename):
    returnVect = zeros((1, 1024))
    fr = open(filename)
    for i in range(32):
        lineStr = fr.readline()
        for j in range(32):
            returnVect[0, 32*i + j] = int(lineStr[j])
    return returnVect

第七章: 利用AdaBoost元算法提高分类性能

优点: 泛化错误率低,易编码,可以应用再大部分分类器上,无参数调整。
缺点: 对离群点敏感。
使用数据类型: 数值型和标称型数据。

AdaBoost的一般流程

  1. 收集数据: 可以使用任何方法。
  2. 准备数据: 依赖于所使用的弱分类类型,本章使用的是单层决策树,这种分类器可以处理任何数据类型。当然也可以使用任意分类器作为弱分类器,第二章到第六章中的任一分类器都可以充当弱分类器。作为弱分类器,简单分类器的效果更好。
  3. 分析数据: 可以使用任意方法。
  4. 训练算法: AdaBoost的大部分时间都在训练上,分类器将多次在同一数据集上训练弱分类器。
  5. 测试算法: 计算分类的错误率。
  6. 使用算法: 同SVM,AdaBoost预测两个类别中的一个。

示例:在一个数据集上应用AdaBoost

  1. 收集数据: 提供的文本文件
  2. 准备数据: 确保类别标签是+1和-1而非1和0.
  3. 分析数据: 手工检查数据。
  4. 训练算法: 在数据上,利用adaBoostTrainDS()函数训练出一系列的分类器。
  5. 测试算法: 在不采取随机抽样的方法下,对AdaBoost和Logistic回归的结果进行完全对等的比较。
  6. 使用算法: 观察例子上的错误率。

adaboost.py完整代码(章节中Python命令行代码在GitHub):

from numpy import *
import matplotlib.pyplot as plt

def loadSimpData():
    datMat = matrix([[1., 2.1],
                        [2. , 1.1],
                        [1.3, 1. ],
                        [1. , 1. ],
                        [2. , 1. ]])
    classLabels = [1.0, 1.0, -1.0, -1.0, 1.0]
    return datMat, classLabels

# Single layer decision tree (decision stump), program 7_1
def stumpClassify(dataMatrix, dimen, threshVal, threshIneq):
    retArray = ones((shape(dataMatrix)[0], 1))
    if threshIneq == 'lt':
        retArray[dataMatrix[:, dimen] <= threshVal] = -1.0
    else:
        retArray[dataMatrix[:, dimen] > threshVal] = -1.0
    return retArray

def buildStump(dataArr, classLabels, D):
    dataMatrix = mat(dataArr)
    labelMat = mat(classLabels).T
    m, n = shape(dataMatrix)
    numSteps = 10.0
    bestStump = {}
    bestClasEst = mat(zeros((m, 1)))
    minError = inf
    for i in range(n):
        rangeMin = dataMatrix[:, i].min()
        rangeMax = dataMatrix[:, i].max()
        stepSize = (rangeMax - rangeMin) / numSteps
        for j in range(-1, int(numSteps) + 1):
            for inequal in ['lt', 'gt']:
                threshVal = (rangeMin + float(j) * stepSize)
                predictedVals = \
                    stumpClassify(dataMatrix, i, threshVal, inequal)
                errArr = mat(ones((m, 1)))
                errArr[predictedVals == labelMat] = 0
                weightedError = D.T*errArr
                # calculate the weighted error rate
                # print("split: dim %d, thresh %.2f, thresh inequal: \
                #       %s, the weighted error is %.3f" %\
                #       (i, threshVal, inequal, weightedError))
                if weightedError < minError:
                    minError = weightedError
                    bestClasEst = predictedVals.copy()
                    bestStump['dim'] = i
                    bestStump['thresh'] = threshVal
                    bestStump['ineq'] = inequal
    return bestStump, minError, bestClasEst

# AdaBoost training function, base on decision stump, program 7_2
def adaBoostTrainDS(dataArr, classLabels, numIt = 40):
    weakClassArr = []
    m = shape(dataArr)[0]
    D = mat(ones((m, 1)) / m)
    aggClassEst = mat(zeros((m, 1)))
    for i in range(numIt):
        bestStump, error, classEst = buildStump(dataArr, classLabels, D)
        print("D:", D.T)
        alpha = float(0.5 * log((1.0 - error)/ max(error, 1e-16)))
        bestStump['alpha'] = alpha
        weakClassArr.append(bestStump)
        print("classEst: ", classEst.T)
        # Calculate D value for the next iteration
        expon = multiply(-1 * alpha * mat(classLabels).T, classEst)
        D = multiply(D, exp(expon))
        D = D / D.sum()
        # Accumulate the error rate
        aggClassEst += alpha * classEst
        print("aggClassEst: ", aggClassEst.T)
        aggErrors = multiply(sign(aggClassEst) != mat(classLabels).T, \
                                                    ones((m, 1)))
        errorRate = aggErrors.sum() / m
        print("total error: ", errorRate, "\n")
        if errorRate == 0.0:
            break
    # program 7_2 version
    # return weakClassArr
    # program 7_5 version
    return weakClassArr, aggClassEst

# AdaBoost's classification function, program 7_3
def adaClassify(datToClass, classifierArr):
    dataMatrix = mat(datToClass)
    m = shape(dataMatrix)[0]
    aggClassEst = mat(zeros((m, 1)))
    for i in range(len(classifierArr)):
        classEst = stumpClassify(dataMatrix, classifierArr[i]['dim'], \
                                classifierArr[i]['thresh'], \
                                classifierArr[i]['ineq'])
        aggClassEst += classifierArr[i]['alpha'] * classEst
        print(aggClassEst)
    return sign(aggClassEst)

# Self adaptive data loading funtion, program 7_4
def loadDataSet(fileName):
    numFeat = len(open(fileName).readline().split('\t'))
    dataMat = []
    labelMat = []
    fr = open(fileName)
    for line in fr.readlines():
        lineArr = []
        curLine = line.strip().split('\t')
        for i in range(numFeat - 1):
            lineArr.append(float(curLine[i]))
        dataMat.append(lineArr)
        labelMat.append(float(curLine[-1]))
    return dataMat, labelMat

def plotROC(preStrengths, classLabels):
    cur = (1.0, 1.0)
    ySum = 0.0
    numPosClas = sum(array(classLabels) == 1.0)
    yStep = 1 / float(numPosClas)
    xStep = 1 / float(len(classLabels) - numPosClas)
    # Get a sorted list of index
    sortedIndicies = preStrengths.argsort()
    fig = plt.figure()
    fig.clf()
    ax = plt.subplot(111)
    for index in sortedIndicies.tolist()[0]:
        if classLabels[index] == 1.0:
            delX = 0
            delY = yStep
        else:
            delX = xStep
            delY = 0
            ySum += cur[1]
        ax.plot([cur[0], cur[0] - delX], [cur[1], cur[1] - delY], c = 'b')
        cur = (cur[0] - delX, cur[1] - delY)
    ax.plot([0, 1], [0, 1], 'b--')
    plt.xlabel('False Positive Rate')
    plt.ylabel('True Positive Rate')
    ax.axis([0, 1, 0, 1])
    plt.show()
    print("the Area Under the Curve is: ", ySum * xStep)

]]>
<![CDATA[《机器学习实战》笔记 (第三和第四章) (附Python3版代码)(Machine Learning in Action)]]>https://codeproducers.com/20180318/5e0456874b43c3452b5136dcSun, 18 Mar 2018 12:47:22 GMT

前言:

《机器学习实战》笔记 (第三和第四章) (附Python3版代码)(Machine Learning in Action)

《机器学习实战》(Machine Learning in Action)是一本常见的机器学习入门书,书中代码由Python2写成。由于现时Python2已逐渐退出舞台,所以这篇文章将该书的所有代码部分用Python3重写。
代码上传GitHub: https://github.com/kungbob/Machine_Learning_In_Action
原版Python2代码:https://www.manning.com/books/machine-learning-in-action

第三章:决策树

优点: 计算复杂度不高,输出结果易于理解,对中间值的确实不敏感,可以处理不相关特征数据。
缺点:可能会产生过度匹配问题。
适用数据类型:数据型和标称型。

决策树算法原理:

  1. 寻找当前数据集中的最好特征
  2. 使用该特征划分数据集
  3. 为数据集创建分支节点
  4. 如果数据子集的数据属于同一类型,分类完成。相反,则需要为子集重复整个算法流程。

决策树的一般流程:

  1. 收集数据:可以使用任何方法
  2. 准备数据:树构造算法只适用于标称型数据,因此数值型数据必须离散化。
  3. 分析数据:可以使用任何方法,构造书完成后,应该检查图形是否符合预期。
  4. 训练算法:构造树的数据结构。
  5. 测试算法:使用经验树计算错误率。
  6. 使用算法:此步骤可以适用于任何监督学习算法,而使用决策树可以更好地理解数据的内在含义。

注: 书中使用的是ID3算法来划分数据,并非二分法。

示例:使用决策树预测隐形眼镜类型

  1. 收集数据:提供的文本文件。
  2. 准备数据:解析tab键分隔的数据行。
  3. 分析数据:快速检查数据,确保正确地解析数据内容,使用createPlot()函数绘制最终的树形图。
  4. 训练算法:使用3.1节的createTree()函数。
  5. 测试算法:编写测试函数验证决策树可以正确分类给定的数据实例。
  6. 使用算法:存储树的数据结构,以便下次使用时无需重新构造书。

trees.py完整代码(章节中Python命令行代码在GitHub):

from math import log
import operator

# Entropy calculation, program 3_1
def calcShannonEnt(dataSet):
    numEntries = len(dataSet)
    labelCounts = {}
    # Create dictionary for all possible classes
    for featVec in dataSet:
        currentLabel = featVec[-1]
        if currentLabel not in labelCounts.keys():
            labelCounts[currentLabel] = 0
        labelCounts[currentLabel] += 1
    shannonEnt = 0.0
    # Choose 2 as log number
    for key in labelCounts:
        prob = float(labelCounts[key]) / numEntries
        shannonEnt -= prob * log(prob, 2)
    return shannonEnt

def createDataSet():
    dataSet = [[1, 1, 'yes'],
            [1, 1, 'yes'],
            [1, 0, 'no'],
            [0, 1, 'no'],
            [0, 1, 'no']]
    labels = ['no surfacing', 'flippers']
    return dataSet, labels

# Spliting classes for particular feature, program 3_2
def splitDataSet(dataSet, axis, value):
    # Create a new list object
    retDataSet = []
    for featVec in dataSet:
        if featVec[axis] == value:
            # Extraction of data
            reducedFeatVec = featVec[:axis]
            reducedFeatVec.extend(featVec[axis+1:])
            retDataSet.append(reducedFeatVec)
    return retDataSet

# Finding the best feature to split, program 3_3
def chooseBestFeatureToSplit(dataSet):
    numFeatures = len(dataSet[0]) - 1
    baseEntopy = calcShannonEnt(dataSet)
    bestInfoGain = 0.0
    bestFeature = -1
    for i in range(numFeatures):
        # Create the only one classifier label
        featList = [example[i] for example in dataSet]
        uniqueVals = set(featList)
        newEntropy = 0.0
        # Calculate the information entropy for each classification
        for value in uniqueVals:
            subDataSet = splitDataSet(dataSet, i, value)
            prob = len(subDataSet) / float(len(dataSet))
            newEntropy += prob * calcShannonEnt(subDataSet)
        infoGain = baseEntopy - newEntropy
        # Choose the only best entropy
        if(infoGain > bestInfoGain):
            bestInfoGain = infoGain
            bestFeature = i
    return bestFeature

# Find the majority prediction of classes
def majorityCnt(classList):
    classCount = {}
    for vote in classList:
        if vote not in classCount.keys():
            classCount[vote] = 0
        classCount[vote] += 1
    sortedClassCount = sorted(classCount.iteritems(), \
        key = operator.itergetter(1), reverse = True)
    return sortedClassCount[0][0]

# Create Trees, program 3_4
def createTree(dataSet, labels):
    classList = [example[-1] for example in dataSet]
    # If the class are the same, stop classifying
    if classList.count(classList[0]) == len(classList):
        return classList[0]
    # Return the highest frequency class after iterating all features
    if len(dataSet[0]) == 1:
        return majorityCnt(classList)
    bestFeat = chooseBestFeatureToSplit(dataSet)
    bestFeatLabel = labels[bestFeat]
    myTree = {bestFeatLabel: {}}
    del(labels[bestFeat])
    # Get the list of all features value
    featValues = [example[bestFeat] for example in dataSet]
    uniqueVals = set(featValues)
    for value in uniqueVals:
        subLabels = labels[:]
        myTree[bestFeatLabel][value] = createTree(splitDataSet \
                            (dataSet, bestFeat, value), subLabels)
    return myTree

# classification function of decision tree, program 3_8
def classify(inputTree, featLabels, testVec):
    firstStr = list(inputTree.keys())[0]
    secondDict = inputTree[firstStr]
    featIndex = featLabels.index(firstStr)
    # Change the labeled string into index
    for key in secondDict.keys():
        if testVec[featIndex] == key:
            if type(secondDict[key]).__name__ == 'dict' :
                classLabel = classify(secondDict[key], featLabels, testVec)
            else:
                classLabel = secondDict[key]
    return classLabel

# Use pickle to store decision tree, program 3_9
def storeTree(inputTree, filename):
    fw = open(filename, 'w')
    pickle.dump(inputTree, fw)
    fw.close()

def grabTree(filename):
    fr = open(filename)
    return pickle.load(fr)

treePlotter.py

import matplotlib.pyplot as plt
import pickle

# Define the figure size and arrow format, program 3_5
decisionNode = dict(boxstyle="sawtooth", fc = "0.8")
leafNode = dict(boxstyle="round4", fc = "0.8")
arrow_args = dict(arrowstyle = "<-")

# Draw comment with arrows
def plotNode(nodeTxt, centerPt, parentPt, nodeType):
    createPlot.ax1.annotate(nodeTxt, xy = parentPt, \
    xycoords = 'axes fraction', xytext = centerPt, \
    textcoords = 'axes fraction', \
    va = "center", ha = "center", bbox = nodeType, \
    arrowprops = arrow_args)

def createPlot():
    fig = plt.figure(1, facecolor = 'white')
    fig.clf()
    createPlot.ax1 = plt.subplot(111, frameon = False)
    plotNode('Decision Node', (0.5, 0.1), (0.1, 0.5), decisionNode)
    plotNode('Branch Node', (0.8, 0.1), (0.3, 0.8), leafNode)
    plt.show()

# Number of leafs of tree, program 3_6
def getNumLeafs(myTree):
    numLeafs = 0
    firstStr = list(myTree.keys())[0]
    secondDict = myTree[firstStr]
    # Testing whether the data type in node is dict
    for key in secondDict.keys():
        if type(secondDict[key]).__name__ == 'dict':
            numLeafs += getNumLeafs(secondDict[key])
        else:
            numLeafs += 1
    return numLeafs

def getTreeDepth(myTree):
    maxDepth = 0
    firstStr = list(myTree.keys())[0]
    secondDict = myTree[firstStr]
    for key in secondDict.keys():
        if type(secondDict[key]).__name__ == 'dict':
            thisDepth = 1 + getTreeDepth(secondDict[key])
        else:
            thisDepth = 1
        if thisDepth > maxDepth:
            maxDepth = thisDepth
    return maxDepth

def retrieveTree(i):
    listOfTrees = [{'no surfacing': {0: 'no', 1: {'flippers': \
                    {0: 'no', 1: 'yes'}}}},
                   {'no surfacing': {0: 'no', 1: {'flippers': \
                    {0: {'head': {0: 'no', 1: 'yes'}}, 1: 'no'}}}}
                   ]
    return listOfTrees[i]

# Fill text message in father nodes, program 3_7
def plotMidText(cntrPt, parentPt, txtString):
    xMid = (parentPt[0] - cntrPt[0]) / 2.0 + cntrPt[0]
    yMid = (parentPt[1] - cntrPt[1]) / 2.0 + cntrPt[1]
    createPlot.ax1.text(xMid, yMid, txtString)

# Calculate the height and width
def plotTree(myTree, parentPt, nodeTxt):
    numLeafs = getNumLeafs(myTree)
    depth = getTreeDepth(myTree)
    firstStr = list(myTree.keys())[0]
    cntrPt = (plotTree.xOff + (1.0 + float(numLeafs)) / 2.0 / plotTree.totalW, \
                plotTree.yOff)
    # Mark the child node text
    plotMidText(cntrPt, parentPt, nodeTxt)
    plotNode(firstStr, cntrPt, parentPt, decisionNode)
    secondDict = myTree[firstStr]
    # Reduce the shift value of y axis
    plotTree.yOff = plotTree.yOff - 1.0/plotTree.totalD
    for key in secondDict.keys():
        if type(secondDict[key]).__name__== 'dict':
            plotTree(secondDict[key], cntrPt, str(key))
        else:
            plotTree.xOff = plotTree.xOff + 1.0 / plotTree.totalW
            plotNode(secondDict[key], (plotTree.xOff, plotTree.yOff), \
                cntrPt, leafNode)
            plotMidText((plotTree.xOff, plotTree.yOff), cntrPt, str(key))
    plotTree.yOff = plotTree.yOff + 1.0 / plotTree.totalD

def createPlot(inTree):
    fig = plt.figure(1, facecolor = 'white')
    fig.clf()
    axprops = dict(xticks = [], yticks = [])
    createPlot.ax1 = plt.subplot(111, frameon = False, **axprops)
    plotTree.totalW = float(getNumLeafs(inTree))
    plotTree.totalD = float(getTreeDepth(inTree))
    plotTree.xOff = -0.5 / plotTree.totalW
    plotTree.yOff = 1.0
    plotTree(inTree, (0.5, 1.0), '')
    plt.show()

第四章:基于概率论的分类方法:朴素贝叶斯

优点: 在数据较少的情况下仍然有效,可以处理多类别问题。
缺点: 对于输入数据的准备方式较为敏感。
适用数据类型: 标称型数据。

朴素贝叶斯的一般过程:

  1. 收集数据: 可以使用任何方法。本章使用RSS源。
  2. 准备数据: 需要数值型或者布尔型数据。
  3. 分析数据: 有大量特征时,绘制特征作用不大,此时使用时直方图效果更好。
  4. 训练算法: 计算不同的独立特征的条件概率。
  5. 测试算法: 计算错误率。
  6. 使用算法: 一个常见的朴素贝叶斯应用是文档分类。可以在任意的分类场景中使用朴素贝叶斯分类器,不一定非要是文本。

示例:使用朴素贝叶斯对电子邮件进行分类

  1. 收集数据: 提供文本文件。
  2. 准备数据: 将文本文件解析成词条向量。
  3. 分析数据: 检查词条确保解析的正确性。
  4. 训练算法: 使用之前建立的trainNB0()函数。
  5. 测试算法: 使用classifyNB(),并且构建一个新的测试函数来计算文档集的错误率。
  6. 使用算法: 构建一个完整的程序对一组文档进行分类,将错分的文档输出到屏幕上。

bayes.py完整代码(章节中Python命令行代码在GitHub):

from numpy import *
import re, operator

# Convert wordlist to vector, program 4_1
def loadDataSet():
    postingList=[['my', 'dog', 'has', 'flea', 'problems', 'help', 'please'],
                 ['maybe', 'not', 'take', 'him', 'to', 'dog', 'park', 'stupid'],
                 ['my', 'dalmation', 'is', 'so', 'cute', 'I', 'love', 'him'],
                 ['stop', 'posting', 'stupid', 'worthless', 'garbage'],
                 ['mr', 'licks', 'ate', 'my', 'steak', 'how', 'to', 'stop', 'him'],
                 ['quit', 'buying', 'worthless', 'dog', 'food', 'stupid']]
    classVec = [0,1,0,1,0,1]    #1 is abusive, 0 not
    return postingList,classVec

def createVocabList(dataSet):
    # Create an empty set
    vocabSet = set([])
    for document in dataSet:
        # Create a union set
        vocabSet = vocabSet | set(document)
    return list(vocabSet)


def setOfWords2Vec(vocabList, inputSet):
    # Create a zero vector
    returnVec = [0] * len(vocabList)
    for word in inputSet:
        if word in vocabList:
            returnVec[vocabList.index(word)] = 1
        else:
            print("the word: %s is not in mu vocabulary!" % word)
    return returnVec

# Naive Bayes classification training function, program 4_2
def trainNB0(trainMatrix, trainCategory):
    numTrainDocs = len(trainMatrix)
    numWords = len(trainMatrix[0])
    pAbusive = sum(trainCategory)/float(numTrainDocs)
    # Initialise the probability as zero in program 4_2
    # Change in chapter 4_5_3
    p0Num = ones(numWords)
    p1Num = ones(numWords)
    p0Denom = 2.0
    p1Denom = 2.0
    for i in range(numTrainDocs):
        if trainCategory[i] == 1:
            # Addition of vectors
            p1Num += trainMatrix[i]
            p1Denom += sum(trainMatrix[i])
        else:
            p0Num += trainMatrix[i]
            p0Denom += sum(trainMatrix[i])
    # Division to each element only in program 4_2
    # Add log function in chapter 4_5_3
    p1Vect = log(p1Num/p1Denom) # change to log()
    p0Vect = log(p0Num/p0Denom) # change to log()
    return p0Vect, p1Vect, pAbusive

# Naive Bayes classification function
# Program 4_3
def classifyNB(vec2Classify, p0Vec, p1Vec, pClass1):
    p1 = sum(vec2Classify * p1Vec) + log(pClass1)
    p0 = sum(vec2Classify * p0Vec) + log(1.0 - pClass1)
    if p1 > p0:
        return 1
    else:
        return 0

def testingNB():
    listOPosts, listClasses = loadDataSet()
    myVocabList = createVocabList(listOPosts)
    trainMat = []
    for postinDoc in listOPosts:
        trainMat.append(setOfWords2Vec(myVocabList, postinDoc))
    p0V, p1V, pAb = trainNB0(array(trainMat), array(listClasses))
    testEntry = ['love', 'my', 'dalmation']
    thisDoc = array(setOfWords2Vec(myVocabList, testEntry))
    print(testEntry, 'classified as: ', classifyNB(thisDoc, p0V, p1V, pAb))
    testEntry = ['stupid', 'garbage']
    thisDoc = array(setOfWords2Vec(myVocabList, testEntry))
    print(testEntry, 'classified as: ', classifyNB(thisDoc, p0V, p1V, pAb))

# Naive Bayes Word Bag Model, program 4_4
def bagOfWords2VecMN(vocabList, inputSet):
    returnVec = [0] * len(vocabList)
    for word in inputSet:
        if word in vocabList:
            returnVec[vocabList.index(word)] += 1
    return returnVec

# File analysis and complte email test function, program 4_5
def textParse(bigString):
    listOfTokens = re.split(r'\W*', bigString)
    return [tok.lower() for tok in listOfTokens if len(tok) > 2]

def spamTest():
    docList = [];
    classList = [];
    fullText = [];
    for i in range(1, 26):
        # Load and parse the email
        wordList = textParse(open('email/spam/%d.txt' %i).read())
        docList.append(wordList)
        fullText.extend(wordList)
        classList.append(1)
        wordList = textParse(open('email/ham/%d.txt' %i).read())
        docList.append(wordList)
        fullText.extend(wordList)
        classList.append(0)
    vocabList = createVocabList(docList)
    trainingSet = list(range(50))
    testSet = []
    for i in range(10):
        # Construct training dataset using random number
        randIndex = int(random.uniform(0, len(trainingSet)))
        testSet.append(trainingSet[randIndex])
        del(trainingSet[randIndex])
    trainMat = []
    trainClasses = []
    for docIndex in trainingSet:
        trainMat.append(setOfWords2Vec(vocabList, docList[docIndex]))
        trainClasses.append(classList[docIndex])
    p0V, p1V, pSpam = trainNB0(array(trainMat), array(trainClasses))
    errorCount = 0
    for docIndex in testSet:
        # Classify the training dataset
        wordVector = setOfWords2Vec(vocabList, docList[docIndex])
        if classifyNB(array(wordVector), p0V, p1V, pSpam) != classList[docIndex]:
            errorCount += 1
    print('the error rate is: ', float(errorCount)/len(testSet))

# RSS classifier and high drequency word deletion, program 4_6
def calcMostFreq(vocabList, fullText):
    # Calculate the frequency of word
    freqDict = {}
    for token in vocabList:
        freqDict[token] = fullText.count(token)
    sortedFreq = sorted(freqDict.items(), key=operator.itemgetter(1), \
                    reverse = True)
    return sortedFreq[:30]

def localWords(feed1, feed0):
    docList = []
    classList = []
    fullText = []
    minLen = min(len(feed1['entries']), len(feed0['entries']))
    for i in range(minLen):
        # Visit one RSS source each time
        wordList = textParse(feed1['entries'][i]['summary'])
        docList.append(wordList)
        fullText.extend(wordList)
        classList.append(1)
        wordList = textParse(feed0['entries'][i]['summary'])
        docList.append(wordList)
        fullText.extend(wordList)
        classList.append(0)
    vocabList = createVocabList(docList)
    # Remove the top 30 words from vocabList
    top30Words = calcMostFreq(vocabList, fullText)
    for pairW in top30Words:
        if pairW[0] in vocabList:
            vocabList.remove(pairW[0])
    trainingSet = list(range(2*minLen))
    testSet = []
    for i in range(20):
        randIndex = int(random.uniform(0, len(trainingSet)))
        testSet.append(trainingSet[randIndex])
        del(trainingSet[randIndex])
    trainMat = []
    trainClasses = []
    for docIndex in trainingSet:
        trainMat.append(bagOfWords2VecMN(vocabList, docList[docIndex]))
        trainClasses.append(classList[docIndex])
    p0V, p1V, pSpam = trainNB0(array(trainMat), array(trainClasses))
    errorCount = 0
    for docIndex in testSet:
        wordVector = bagOfWords2VecMN(vocabList, docList[docIndex])
        if classifyNB(array(wordVector), p0V, p1V, pSpam) != classList[docIndex]:
            errorCount += 1
    print ('the error rate is: ', float(errorCount) / len(testSet))
    return vocabList, p0V, p1V

# Most frequency word showing function
def getTopWords(ny, sf):
    vocabList, p0V, p1V = localWords(ny, sf)
    topNY = []
    topSF = []
    for i in range(len(p0V)):
        if p0V[i] > -6.0:
            topSF.append((vocabList[i], p0V[i]))
        if p1V[i] > -6.0:
            topNY.append((vocabList[i], p1V[i]))
    sortedSF = sorted(topSF, key=lambda pair: pair[1], reverse = True)
    print("SF**SF**SF**SF**SF**SF**SF**SF**SF**SF**SF**SF**SF**")
    for item in sortedSF:
        print(item[0])
    sortedNY = sorted(topNY, key=lambda pair: pair[1], reverse = True)
    print("NY**NY**NY**NY**NY**NY**NY**NY**NY**NY**NY**NY**NY**")
    for item in sortedNY:
        print(item[0])

示例: 使用朴素贝叶斯分类器从个人广告中获取区域倾向

注意: 现时书中的RSS网址已经不可用,暂时这部分代码并不能运行。

  1. 收集数据: 从RSS源收集内容,这里需要对RSS源构建一个接口。
  2. 准备数据: 将文本文件解构成词条向量。
  3. 分析数据: 检查词条确保解析的正确性。
  4. 训练算法: 使用之前建立的trainNB0()函数。
  5. 测试算法: 观察错误率,确保分类器可用。可以修改切分程序,以降低错误率,提高分类结果。
  6. 使用算法: 构建一个完整的程序,封装所有内容。给定两个RSS源,该程序会显示最常用的公共词。

该部分的代码已经整合于bayes.py文件中。

]]>
<![CDATA[《机器学习实战》笔记 (第一和第二章) (附Python3版代码)(Machine Learning in Action)]]>https://codeproducers.com/20180115/5e0456874b43c3452b5136dbMon, 15 Jan 2018 13:46:08 GMT

前言:

《机器学习实战》笔记 (第一和第二章) (附Python3版代码)(Machine Learning in Action)

《机器学习实战》(Machine Learning in Action)是一本常见的机器学习入门书,书中代码由Python2写成。由于现时Python2已逐渐退出舞台,所以这篇文章将该书的所有代码部分用Python3重写。
代码上传GitHub: https://github.com/kungbob/Machine_Learning_In_Action
原版Python2代码:https://www.manning.com/books/machine-learning-in-action

第一章:机器学习基础

机器学习一个主要任务是分类,另外一项是回归(预测数据)。

开发机器学习应用程序的步骤:

  1. 收集数据
  2. 准备输入数据
  3. 分析输入数据
  4. 训练算法
  5. 测试算法
  6. 使用算法
# import numpy library
from numpy import *
# randomly generate four lists of four elements
print("Matrix is:\n", random.rand(4, 4))
# change the list into a 4 by 4 matrix
randMat = mat(random.rand(4, 4))
# calculate the reverse of the matrix
invRandMat = randMat.I
print("Inverse is:\n", invRandMat.I)
# Matrix * Inverse = Identity Matrix
myEye = randMat * invRandMat
print("Matrix * Inverse is:\n", myEye)
# slight difference with real Identity Matrix deal to computation error
diff = myEye - eye(4)
print("Difference is:\n", diff)

第二章:k-近邻算法

优点: 精度高,对异常值不敏感,无数据输入假定。
缺点:计算复杂度高,空间复杂负高。
适用数据范围:数值型和标称性。

k-近邻算法(kNN)算法原理:

  1. 根据训练样本集的标签进行分类,知道样本和所属分类的对应关系。
  2. 输入没有标签的新数据,将新数据的每个特征与样本集数据的对应特征进行比较。
  3. 实用算法提取样本中最相似数据的分类标签。
  4. k代表的就是数据集中k个最相似的数据(通常k不大于20)。
  5. 根据k个最相似数据中出现次数最多的分类,作为新数据的分类。

k-近邻算法的一般流程:

  1. 收集数据:任何方法。
  2. 准备数据:距离计算所需要的数值,最好是结构化的数据格式。
  3. 分析数据:任何方法。
  4. 训练算法:此步骤不适用于k-近邻算法。
  5. 测试算法:计算错误率。
  6. 使用算法:首先需要输入样本数据和结构化的结果,然后运行k-近邻算法判定输入数据分别属于哪个类别,最后应用对计算出的分类执行后续的处理。

示例:手写识别系统

  1. 收集数据:提供文本文件。
  2. 准备数据:编写函数img2vector(),将图像格式转换为分类器使用的向量格式。
  3. 分析数据:在Python命令提示符中检查数据,确保符合要求。
  4. 训练算法:此步骤不适用于k-近邻算法。
  5. 测试算法:编写函数使用提供的部分数据集作为测试样本,测试样本与非测试样本的却别在于测试样本是已经完成分类的数据,如果预测分类与实际分类不同,则标记为一个错误。
  6. 使用算法:书中没有完成此步骤,读者感兴趣可以构建完整的应用程序,从图像中提取数字,完成数字识别。

kNN.py完整代码(章节中Python命令行代码在GitHub):

from numpy import *
from os import listdir
import operator

# Create data for later usage, chapter 2_1_1
def createDataSet():
    group = array([[1.0, 1.1], [1.0, 1.0], [0, 0], [0, 0.1]])
    labels = ['A', 'A', 'B', 'B']
    return group, labels

# kNN algorithm, program 2_1
def classify0(inX, dataSet, labels, k):
    dataSetSize = dataSet.shape[0]
    # Distance between two points, using eculidean
    diffMat = tile(inX, (dataSetSize, 1)) - dataSet
    sqDiffMat = diffMat**2
    sqDistances = sqDiffMat.sum(axis = 1)
    distances = sqDistances**0.5
    # Finding the k point of shortest distance
    sortedDistIndices = distances.argsort()
    classCount = {}
    for i in range(k):
        voteIlabel = labels[sortedDistIndices[i]]
        classCount[voteIlabel] = classCount.get(voteIlabel, 0) + 1
    # Sorting
    sortedClassCount = sorted(classCount.items(),
      key = operator.itemgetter(1), reverse = True)
    return sortedClassCount[0][0]

# Read txt file into matrix, program 2_2
def file2matrix(filename):
    fr = open(filename)
    arrayOLines = fr.readlines()
    # Get the lines of file
    numberOfLines = len(arrayOLines)
    # Create a matrix with all zeros
    returnMat = zeros((numberOfLines, 3))
    classLabelVector = []
    index = 0
    # Read numbers into matrix
    for line in arrayOLines:
        line = line.strip()
        listFromLine = line.split('\t')
        returnMat[index, :] = listFromLine[0:3]
        classLabelVector.append(int(listFromLine[-1]))
        index += 1
    return returnMat, classLabelVector

# Normalization of features, program 2_3
def autoNorm(dataSet):
    minVals = dataSet.min(0)
    maxVals = dataSet.max(0)
    ranges = maxVals - minVals
    normDataSet = zeros(shape(dataSet))
    m = dataSet.shape[0]
    # Normalization Step, dataset minus minimum value
    # tile function means create a matrix with m row, repeating minVals for
    # column once.
    normDataSet = dataSet - tile(minVals, (m, 1))
    # Normalization Step, divide the range of value
    normDataSet = normDataSet / tile(ranges, (m, 1))
    return normDataSet, ranges, minVals

# Classifer for dating website, testing its accuracy, program 2_4
def datingClassTest():
    # Ratio of testing case
    hoRatio = 0.1
    datingDataMat, datingLabels = file2matrix('datingTestSet2.txt')
    normMat, ranges, minVals = autoNorm(datingDataMat)
    m = normMat.shape[0]
    numTestVecs = int(m * hoRatio)
    errorCount = 0.0
    # Testing classifier
    for i in range(numTestVecs):
        classifierResult = classify0(normMat[i, :], normMat[numTestVecs:m, :],
                            datingLabels[numTestVecs:m], 3)
        print("The classifier came back with: %d, the real answer is: %d" \
                % (classifierResult, datingLabels[i]))
        if classifierResult != datingLabels[i] :
            errorCount += 1.0
    print("The total error rate is: %f" %(errorCount/float(numTestVecs)))

# Dating Website Prediction Function, program 2_5
def classifyPerson():
    resultList = ['not at all', 'in small doses', 'in large doses']
    percentTats = float(input(
                    "percentage of time spent playing video games?"))
    ffMiles = float(input("frequent flier miles earned per year?"))
    iceCream = float(input("liters of ice cream consumed per year?"))
    datingDataMat, datingLabels = file2matrix('datingTestSet2.txt')
    normMat, ranges, minVals = autoNorm(datingDataMat)
    inArr = array([ffMiles, percentTats, iceCream])
    classifierResult = classify0((inArr - \
                        minVals) / ranges, normMat, datingLabels, 3)
    print("You will probably like this person: ", \
                    resultList[classifierResult - 1])

# Change image of 32 x 32 pixels to 1 x 1024 vector, chapter 2_3_1
def img2vector(filename):
    returnVect = zeros((1, 1024))
    fr = open(filename)
    for i in range(32):
        lineStr = fr.readline()
        for j in range(32):
            returnVect[0, 32*i + j] = int(lineStr[j])
    return returnVect

# Handwriting Testing Function, chapter 2_6
def handwritingClassTest():
    # Retrieve directory content
    hwLabels = []
    trainingFileList = listdir('trainingDigits')
    m = len(trainingFileList)
    trainingMat = zeros((m, 1024))
    # Analysis the number from the file name, training the data
    for i in range(m):
        fileNameStr = trainingFileList[i]
        fileStr = fileNameStr.split('.')[0]
        classNumStr = int(fileStr.split('_')[0])
        hwLabels.append(classNumStr)
        trainingMat[i, :] = img2vector('trainingDigits/%s' % fileNameStr)
    testFileList = listdir('testDigits')
    errorCount = 0.0
    mTest = len(testFileList)
    # Testing the algorithm
    for i in range(mTest):
        fileNameStr = testFileList[i]
        fileStr = fileNameStr.split('_')[0]
        classNumStr = int(fileStr.split('_')[0])
        vectorUnderTest = img2vector('testDigits/%s' % fileNameStr)
        classifierResult = classify0(vectorUnderTest, \
                                trainingMat, hwLabels, 3)
        print("The classifier came back with: %d, the real answer is: %d"\
                % (classifierResult, classNumStr))
        if classifierResult != classNumStr:
            errorCount += 1.0
    print ("\nthe total number of errors is: %d" % errorCount)
    print ("\nthe total error rate is: %f" % (errorCount/float(mTest)))
]]>
<![CDATA[Andrew Ng机器学习课程笔记(第八周至第十一周)]]>https://codeproducers.com/20171229-2/5e0456874b43c3452b5136daFri, 29 Dec 2017 08:00:16 GMTAndrew Ng机器学习课程笔记(第八周至第十一周)

本笔记基于Coursera上吴恩达(Andrew Ng)的机器学习课程内容集结而成。(课程链接)
全课程分为十一周,每周都有1-3个测验。第二周至第九周附有编程作业,一共8份。

更新:
3/1/2018: 更正kNN为k-平均算法

第八周

非监督学习(Unsupervised Learning)

非监督学习是使用非标签的数据来建立模型。非监督学习可用于建立数据集群,为数据分类。

k-平均演算法 (K-Means Algoithm)

  1. 随机在数据集中选择两点为集群中心(cluster centroids)。
  2. 为所有数据选择最近的集群中心。
  3. 移动集群中心到下属数据的中心点。
  4. 重新运行2-3直至建立集群。

K代表的是集群的数量。\({x^{(1)}, x^{(2)}, \dots,x^{(m)}}\)代表训练组。
在为数据选择最近的集群中心的步骤中,公式为:
\[c^{(i)} = argmin_k\ ||x^{(i)} - \mu_k||^2\]
其中\(c^{(i)}\)代表的是离数据\(x^{(i)}\)最近的集群中心的编号。
移动集群中心到下属数据的中心点,公式为:
\[\mu_k = \dfrac{1}{n}[x^{(k_1)} + x^{(k_2)} + \dots + x^{(k_n)}] \in \mathbb{R}^n\]
如果一个集群中心并没有任何下属的数据,可以为该中心重新随机设计一点,也可以直接消除该数据集。在一定数量的循环过后,模型会收敛,集群并不会有所变化。

优化目标

代价函数:
\[J(c^{(i)},\dots,c^{(m)},\mu_1,\dots,\mu_K) = \dfrac{1}{m}\sum_{i=1}^m ||x^{(i)} - \mu_{c^{(i)}}||^2\]
其中,\(c^{(i)}\)代表的是数据\(x^{(i)}\)所属的集群中心编号;\(\mu_{(k)}\)是集群中心;\(\mu_{c^(i)}\)代表的是\(x^{(i)}\)所属的集群中心。
优化目标是要\(min_{c,\mu}\ J(c,\mu)\)。代价函数的数值只会下降,不会上升。

随机初始

  1. 保证集群的数量少于训练样本数量。(K < m)
  2. 随机选择K个不重复的训练样本。
  3. 设置\(\mu_1,\dots,\mu_K\)在训练样本上。

k-平均演算法可能会困在局部最优解。为了减少发生的几率,算法模型需要至少跑几次以上。

集群的数量

肘部法则:为Cost J和K画图显示两者关系。Cost J在K到一定数量后,下降的幅度减少,趋于平缓。找出那一点的K值。(通常是Cost J的曲线是逐步下降的,那一转折点并不明显。)

下游目的(downstream purpose):选择一个最能达到想要目地的K值。也就是说,在选择K个集群,出来的结果最为有用,能达到想要的目的。

数据降维

降到数据的维度,以达到减少臃肿数据的目的。同时,透过该方法也能做到数据可视化的目的。

主成份分析(Principal Component Analysis (PCA))

PCA的目的是要减少每一个特征到投影线的平均距离,也就是减少投影误差。
PCA并不是线性回归,线性回归是减少平方误差,PCA是减少垂直距离。
预先处理:

  1. 给予一个训练集\(x^{(1)},x^{(2)},x^{(3)}...x(m)\)
  2. 特征缩放\(\mu_j = \dfrac{1}{m}\sum^m_{i=1}x_j^{(i)}\)
  3. 用\(x_j^{(i)} - \mu_j\)取代\(x_j^{(i)}\)
  4. 缩放所有特征\(x_j^{(i)} = \dfrac{x_j^{(i)} - \mu_j}{s_j}\)

2D数据缩放到1D数据,可以用公式来表示:
\[\Sigma = \dfrac{1}{m}\sum^m_{i=1}(x^{(i)})(x^{(i)})^T\]

PCA过程:

  1. 计算出协方差矩阵(covariance matrix)\(\Sigma = \dfrac{1}{m}\sum^m_{i=1}(x^{(i)})(x^{(i)})^T\)
  2. 找出协方差矩阵的特征向量(eigenvector)
  3. 取出头K个栏位的矩阵计算出z \(z^{(i)} = Ureduce^T \cdot x^{(i)}\)
    总结所有步骤:
Sigma = (1/m) * X' * X; % compute the covariance matrix
[U,S,V] = svd(Sigma);   % compute our projected directions
Ureduce = U(:,1:k);     % take the first k directions
Z = X * Ureduce;        % compute the projected data points

由压缩数据重建,由1D变2D的话(\(z \in \mathbb{R} \rightarrow x \in \mathbb{R}^2\)):\(x_{approx}^{(1)} = U_{reduce} \cdot z^{(1)}\)

K的数值可以参考这方法:

  1. 计算出平均平方投影错误。\(\dfrac{1}{m}\sum^m_{i=1}||x^{(i)} - x_{approx}^{(i)}||^2\)
  2. 计算出\(\dfrac{1}{m}\sum^m_{i=1}||x^{(i)}||^2\)
  3. 选一个最小的k值
    \[\dfrac{\dfrac{1}{m}\sum^m_{i=1}||x^{(i)} - x_{approx}^{(i)}||^2}{\dfrac{1}{m}\sum^m_{i=1}||x^{(i)}||^2} \leq 0.01\]
    可以这样表示:
    \[\dfrac{\sum_{i=1}^kS_{ii}}{\sum_{i=1}^nS_{ii}} \geq 0.99\]

使用PCA来减少过拟合问题是一个坏的用法。不要假设每个模型都需要PCA,先正常运行一次,再考虑是否需要PCA。

第九周

异常检测

设置\(x^{(i)}\)为i的特征。用数据设置模型p(x)。根据\(p(x) < ϵ\)识别一些不正常的活动。如果检测出过多/太多的异常,可以减少阀值\(ϵ\)。

高斯分布(正态分布)

x的高斯机会分布可以这样表示:\(x \sim \mathcal{N}(\mu, \sigma^2)\)。
完整的公式表示是:
\[\large p(x;\mu,\sigma^2) = \dfrac{1}{\sigma\sqrt{(2\pi)}}e^{-\dfrac{1}{2}(\dfrac{x - \mu}{\sigma})^2}\]
估算参数\(\mu\)为数据集的平均值:
\[\mu = \dfrac{1}{m}\displaystyle \sum_{i=1}^m x^{(i)}\]
然后估算参数\(\sigma^2\),用平均误差公式:
\[\sigma^2 = \dfrac{1}{m}\displaystyle \sum_{i=1}^m(x^{(i)} - \mu)^2\]
假设数据集的每个数据都是独立的,我们可以得出公式:
\[p(x) = \displaystyle \prod^n_{j=1} p(x_j;\mu_j,\sigma_j^2) \]

评估检测系统

用训练集\(\lbrace x^{(1)},\dots,x^{(m)} \rbrace\)算出模型系统的\(p(x)\)值。
在验证集上,预测一下结果:
如果 \(p(x) < ϵ\),y = 1
如果 \(p(x) \ge ϵ\),y = 0
然后使用F1数值等方法评估模型。

使用异常检测的情况:

  • 有一个很小数目的正样本和很多数目的负样本
  • 有很多种不同的异常,难于从正样本中找出,未来异常可能跟以往的异常情况不一样。

使用监督学习的情况:

  • 大数目的政府样本
  • 有足够多的正样本,可以推测出新的正样本的样子。

可以为数据进行变形(transforms)以便适合高斯分布。(例如:\(\log(x),\log(x+1),\sqrt(x),x^{1/3}\))

多元高斯分布
\[p(x;\mu,\Sigma) = \dfrac{1}{(2\pi)^{n/2} |\Sigma|^{1/2}} exp(-1/2(x-\mu)^T\Sigma^{-1}(x-\mu))\]

推荐系统

\(n_u \)为用户的数量,\(n_m \)为电影的数量,\(r(i,j) = 1\)为用户j评价了电影i,\(y(i,j) = 1\)为用户j给电影i的评分。
\(\theta^{(j)} \)为用户j的参数向量,\(x^{(i)} \)为电影i的参数向量。预测用户j在电影i的评分会是:\((\theta^{(j)})^Tx^{(i)}\)。
设置\(m^{(j)}\)为用户j评分过的电影数量。
\(\theta^{(j)}\)的数值可以透过这方法得出:
\[min_{\theta^{(1)},\dots,\theta^{(n_u)}} = \dfrac{1}{2}\displaystyle \sum_{j=1}^{n_u} \sum_{i:r(i,j)=1} ((\theta^{(j)})^T(x^{(i)}) - y^{(i,j)})^2 + \dfrac{\lambda}{2} \sum_{j=1}^{n_u} \sum_{k=1}^n(\theta_k^{(j)})^2\]

协同过滤(collaborative filtering)

代价函数可以更改为:
\[J(x,\theta) = \dfrac{1}{2} \displaystyle \sum_{(i,j):r(i,j)=1}((\theta^{(j)})^Tx^{(i)} - y^{(i,j)})^2 + \dfrac{\lambda}{2}\sum_{i=1}^{n_m} \sum_{k=1}^{n} (x_k^{(i)})^2 + \dfrac{\lambda}{2}\sum_{j=1}^{n_u} \sum_{k=1}^{n} (\theta_k^{(j)})^2\]
算法分为几个步骤:

  1. 初始设置\(x^{(i)},...,x^{(n_m)},\theta^{(1)},...,\theta^{(n_u)}\)为一些较小的随机数值。
  2. 用梯度下降的方法最小化\(J(x^{(i)},...,x^{(n_m)},\theta^{(1)},...,\theta^{(n_u)})\)数值。公式为:
    \[\begin{gather}
    x_k^{(i)} := x_k^{(i)} - \alpha\left (\displaystyle \sum_{j:r(i,j)=1}{((\theta^{(j)})^T x^{(i)} - y^{(i,j)}) \theta_k^{(j)}} + \lambda x_k^{(i)} \right) \cr
    \theta_k^{(j)} := \theta_k^{(j)} - \alpha\left (\displaystyle \sum_{i:r(i,j)=1}{((\theta^{(j)})^T x^{(i)} - y^{(i,j)}) x_k^{(i)}} + \lambda \theta_k^{(j)} \right)
    \end{gather}\]
  3. 用\(\theta\)和特征\(x\),预测出评分\(\theta^Tx\)。

第十周

随机梯度下降法

这方法不用一次性计算出所有数据的梯度,而是运用各个数据来缓慢寻找准确的梯度。
\[cost(\theta,(x^{(i)}, y^{(i)})) = \dfrac{1}{2}(h_{\theta}(x^{(i)}) - y^{(i)})^2\]
训练组的代价函数为:
\[J_{train}(\theta) = \dfrac{1}{m} \displaystyle \sum_{i=1}^m cost(\theta, (x^{(i)}, y^{(i)}))\]

该算法的步骤为:

  1. 随机打乱数据组的排序。
  2. 从第一个数据到最后一个数据,计算\(\Theta_j := \Theta_j - \alpha (h_{\Theta}(x^{(i)}) - y^{(i)}) \cdot x^{(i)}_j\)

小型批量梯度下降

假设每次批量处理10个数据:
\[\theta_j := \theta_j - \alpha \dfrac{1}{10} \displaystyle \sum_{k=i}^{i+9} (h_\theta(x^{(k)}) - y^{(k)})x_j^{(k)}\]

随机梯度下降收敛

较小的学习效率很有可能找到一个更好的结果。所以学习效率代数应该要缓慢下降。公式为:
\[\alpha = \dfrac{const1}{iterationNumber + const2}\]

MapReduce和平行数据处理

MapReduce会讲计算工作拆分,最后合并起来计算。
\[\Theta_j := \Theta_j - \alpha \dfrac{1}{z}(temp_j^{(1)} + temp_j^{(2)} + \cdots + temp_j^{(z)})\]

第十一周

机器学习系统流程:图片识别(Photo OCR)

图片识别流水线:图像->文字检测->文字拆分->文字识别

取得更多数据的方法

  1. 换背景
  2. 扭曲/反转
  3. 添加白噪/噪点

上线分析

分析出改进哪一个部分最具性价比。
假设该部分能100%准确做到想要的效果,最终结果能提升多少。

]]>
<![CDATA[Andrew Ng机器学习课程笔记(第四周至第七周)]]>https://codeproducers.com/20171229/5e0456874b43c3452b5136d9Fri, 29 Dec 2017 02:49:47 GMTAndrew Ng机器学习课程笔记(第四周至第七周)

本笔记基于Coursera上吴恩达(Andrew Ng)的机器学习课程内容集结而成。(课程链接)
全课程分为十一周,每周都有1-3个测验。第二周至第九周附有编程作业,一共8份。

第四周

神经网络

如果假设函数十分复杂,单纯的回归公式将会包含数量庞大的参数。使用神经网络是一个好的解决方法。
逻辑激励函数(sigmoid (logistic) activation function):
\[\frac{1}{1 + e^{-\theta^Tx}}\]
\(a_i^{(j)}\)代表的是第\(j\)层中,\(i\)单位的激励。\(\theta^{(j)} \)代表的是矩阵中控制从第\(j\)层到第\(j+1\)层的投影的比重。
如果神经网络只有一层隐含层,看起来就像这样:
\[\begin{bmatrix}x_0 \newline x_1 \newline x_2 \newline x_3\end{bmatrix}\rightarrow\begin{bmatrix}a_1^{(2)} \newline a_2^{(2)} \newline a_3^{(2)} \newline \end{bmatrix}\rightarrow h_\theta(x) \]
而其中,每个激励节点的数值都是根据公式算出来的:
\[
\begin{align*}
a_1^{(2)} = g(\Theta_{10}^{(1)}x_0 + \Theta_{11}^{(1)}x_1 + \Theta_{12}^{(1)}x_2 + \Theta_{13}^{(1)}x_3) \newline
a_2^{(2)} = g(\Theta_{20}^{(1)}x_0 + \Theta_{21}^{(1)}x_1 + \Theta_{22}^{(1)}x_2 + \Theta_{23}^{(1)}x_3) \newline
a_3^{(2)} = g(\Theta_{30}^{(1)}x_0 + \Theta_{31}^{(1)}x_1 + \Theta_{32}^{(1)}x_2 + \Theta_{33}^{(1)}x_3) \newline
h_\Theta(x) = a_1^{(3)} = g(\Theta_{10}^{(2)}a_0^{(2)} + \Theta_{11}^{(2)}a_1^{(2)} + \Theta_{12}^{(2)}a_2^{(2)} + \Theta_{13}^{(2)}a_3^{(2)}) \newline
\end{align*}
\]
至于每一层的\(\theta^{(j)}\),该矩阵的维度是根据该层和下一层的度数决定的。也就是:
如果一个神经网络的第\(j\)层中有\(s_j\)个单位/节点,第\(j+1\)有\(s_{j+1}\)个单位/节点,参数矩阵\(\theta^{(j)}\)的维度则是\(s_{j+1} \times (s_j + 1) \)。额外+1来自于偏差节点,每一层都要新加一个偏差节点来计算。

引用\(z_k^{(j)}\),
\[
\begin{align*}a_1^{(2)} = g(z_1^{(2)}) \newline a_2^{(2)} = g(z_2^{(2)}) \newline a_3^{(2)} = g(z_3^{(2)}) \newline \end{align*}
\]
换句话说,
\[z_k^{(2)} = \Theta_{k,0}^{(1)}x_0 + \Theta_{k,1}^{(1)}x_1 + \cdots + \Theta_{k,n}^{(1)}x_n\]
用矩阵表示:
\[\begin{align*}x = \begin{bmatrix}x_0 \newline x_1 \newline\cdots \newline x_n\end{bmatrix} &z^{(j)} = \begin{bmatrix}z_1^{(j)} \newline z_2^{(j)} \newline\cdots \newline z_n^{(j)}\end{bmatrix}\end{align*}\]
设置\(x = a^{(1)}\),
\[z^{(j)} = \Theta^{(j-1)}a^{(j-1)}\]
\[h_\Theta(x) = a^{(j+1)} = g(z^{(j+1)})\]

多元分类

神经网络算出来的向量可能是这样:
\[h_\Theta(x) =\begin{bmatrix}0 \newline 0 \newline 1 \newline 0 \newline\end{bmatrix}\]
这时候,结果可以分为这样
\[y^{(i)} =\begin{bmatrix}0 \newline 0 \newline 1 \newline 0 \newline\end{bmatrix},\begin{bmatrix}1 \newline 0 \newline 0 \newline 0 \newline\end{bmatrix},\begin{bmatrix}0 \newline 1 \newline 0 \newline 0 \newline\end{bmatrix}, \begin{bmatrix}0 \newline 0 \newline 0 \newline 1 \newline\end{bmatrix}\]

第五周

神经网络的代价函数

定义L为神经网络总共的层数,\(s_l\)为第l层中的节点数,K为输出的节点数。
神经网络的代价函数为:
Andrew Ng机器学习课程笔记(第四周至第七周)

反向传播算法

反向传播算法的目标是计算出\(min_\theta J(\Theta) \),而就是在什么情况下,代价函数的数值最小。同样,先找出代价函数的偏积分,假设\(\delta_j^{(l)}\)为第l层中节点j的误差。
在最后一层中,可以得知\(\delta^{(L)} = a^{(L)} - y\),L是神经网络中的总层数,\(a^{(L)}\)是激励节点的输出向量。最后该公式只是算出神经网络的结果与实际数值的误差。
每一个节点的误差可以用这个公式来表示:
\[\delta^{(l)} = ((\Theta^{(l)})^T \delta^{(l+1)})\ .*\ g'(z^{(l)})\]

代价函数的偏积分(无视正规化):
\[\dfrac{\partial J(\Theta)}{\partial \Theta_{i,j}^{(l)}} = \frac{1}{m}\sum_{t=1}^m a_j^{(t)(l)} {\delta}_i^{(t)(l+1)}\]

流程:

  1. 设置\(\Delta^{(l)}_{i,j} := 0\)
  2. 从t = 1 到 t = m循环:
    a. 设置\(a^{(1)} := x^{(t)}\)
    b. 计算\(a^{(l)}\)
    c. 用\(y^{(t)}\)计算\(\delta^{(L)} = a^{(L)} - y\)
    d. 利用公式\(\delta^{(l)} = ((\Theta^{(l)})^T \delta^{(l+1)})\ .* a^{(l)}\ .* (1 - a^{(l)})\),计算\(\delta^{(L)}\)
    e. \(\Delta^{(l)} := \Delta^{(l)} + a_j^{(l)} \delta_i^{(l+1)}\)或者使用矩阵\(\Delta^{(l)} := \Delta^{(l)} + \delta^{(l+1)}(a^{(l)})^T\)
    f. \(D^{(l)} := \frac{1}{m} (\Delta^{(l)})\)
    g. \(D^{(l)} := \dfrac{1}{m}\Delta^{(l)}\)

反向传播的代价函数

\[
\begin{gather*}J(\theta) = - \frac{1}{m} \sum_{t=1}^m \sum_{k=1}^K \left[ y^{(t)} \ \log (h_\theta (x^{(t)})) + (1 - y^{(t)})\ \log (1 - h_\theta(x^{(t)}))\right] +
\frac{\lambda}{2m}\sum_{l=1}^{L-1} \sum_{i=1}^{s_l} \sum_{j=1}^{s_l+1} ( \theta_{j,i}^{(l)})^2
\end{gather*}
\]
如果只是单一分类和务实正规化,代价函数可以写成:
\[cost(t) =y^{(t)} \ \log (h_\theta (x^{(t)})) + (1 - y^{(t)})\ \log (1 - h_\theta(x^{(t)}))\]

梯度检查

\[
\dfrac{\partial}{\partial\Theta_j}J(\Theta) \approx \dfrac{J(\Theta_1, \dots, \Theta_j + \epsilon, \dots, \Theta_n) - J(\Theta_1, \dots, \Theta_j - \epsilon, \dots, \Theta_n)}{2\epsilon}
\]

随机初始

模型设置时,\(\Theta\)不能初始为0,需要随机选择一个非零正数。应该设置\(\Theta^{(l)}_{ij}\)的范围为\([-\epsilon,\epsilon]\)。公式为:
\[\begin{gather}
\epsilon = \dfrac{\sqrt{6}}{\sqrt{\mathrm{Loutput} + \mathrm{Linput}}} \cr
\Theta^{(l)} = 2 \epsilon ; \mathrm{rand}(\mathrm{Loutput}, \mathrm{Linput} + 1) - \epsilon
\end{gather}\]

组合起来

  1. 随机选择初始\(\Theta\)的比重。
  2. 用正向传播计算出\(h_{\theta}^{(i)}\)。
  3. 设置代价函数。
  4. 用反向传播算出偏积分。
  5. 用梯度检查的方法确认反向传播无误。
  6. 用梯度下降逐步减少误差和修改\(\Theta\)的比重。

第六周

预测误差可以用以下几种方法修正:

  1. 取得更多训练样本
  2. 使用更小的特征集
  3. 使用更多的特征
  4. 使用多元组合特征
  5. 增加或减少\(\lambda\)

评估假设

将训练样本分类两个组/集,一个是训练组,一个是测试组。
先用训练组来训练模型,然后计算出测试组的误差。

测试误差

线性回归:
\[J_{test}(\Theta) = \dfrac{1}{2m_{test}} \sum_{i=1}^{m_{test}}(h_\Theta(x^{(i)}) - y^{(i)})^2\]
逻辑回归:
\[err(h_\Theta(x),y) =
\begin{matrix}
1 & \mbox{if } h_\Theta(x) \geq 0.5\ and\ y = 0\ or\ h_\Theta(x) < 0.5\ and\ y = 1\newline
0 & \mbox otherwise
\end{matrix}\]
平均测试误差:
\[\text{Test Error} = \dfrac{1}{m_{test}} \sum_{i=1}^{m_{test}} err(h_\Theta(x^{(i)}), y^{(i)})\]

模型选择和训练\验证\测试组

没有验证组(不好的方法,不推荐):

  1. 使用训练组和不同多项式次数(polynomial degree)的来训练模型。
  2. 用模型组来测试,选择最佳的多项式次数。
  3. 用\(J_{test}(\Theta^{(d)})\)找出最后的误差。

使用验证组:
在使用验证组的情况下,数据应分为3组,60%归入训练组,20%归入验证组,20%归入测试组。

  1. 用训练组训练模型。
  2. 使用验证组选择最佳的多项式次数。
  3. 用\(J_{test}(\Theta^{(d)})\)算出最后的误差。

诊断偏差和方差问题

高偏差是欠拟合问题,高方差是过拟合问题。如果选择高多项式次数,测试误差会减少。同时,验证误差会减少到一定程度,然后又上升。

高偏差(欠拟合):\(J_{train}(\Theta)\)和\(J_{CV}(\Theta)\)都很大。同时\(J_{CV}(\Theta) \approx J_{train}(\Theta)\)。
高方差(过拟合):\(J_{train}(\Theta)\)很低,但\(J_{CV}(\Theta)\)会远远高于\(J_{train}(\Theta)\)。

正规化和偏差/方差

正规化参数\(\lambda\)也会造成偏差和方差问题。
\(\lambda\)过大:高偏差(欠拟合),\(J_{train}(\Theta)\)和\(J_{CV}(\Theta)\)都很大。
\(\lambda\)适中:正好,\(J_{CV}(\Theta)\)和\( J_{train}(\Theta)\)都差不多。\(J_{CV}(\Theta) \approx J_{train}(\Theta)\)
\(\lambda\)过小:过方差(过拟合),\(J_{train}(\Theta)\)很低,\(J_{CV}(\Theta)\)很大。
为了选择一个最合适的\(\lambda\)数值,选择做一下步骤:

  1. 创建一个\(\lambda\)列表(从0,0.01,0.01,0.02,0.04...)
  2. 使用各种多项式次数来创建模型。
  3. 使用不同\(\lambda\)来训练模型。
  4. 计算出不同模型,不同\(\lambda\)下的验证组误差。
  5. 选择最低的误差组合。
  6. 用测试组来该组合的误差。

学习曲线

高偏差:
小训练组会造成低/小\(J_{train}(\Theta)\)和高/大\(J_{CV}(\Theta)\)。
大训练组会造成\(J_{train}(\Theta)\)和\(J_{CV}(\Theta)\)都很大,同时\(J_{train}(\Theta) \approx J_{CV}(\Theta)\)。
为高偏差模型提供更多的训练数据更不能有所帮助。

高方差:
小训练组会造成低/小\(J_{train}(\Theta)\)和高/大\(J_{CV}(\Theta)\)。
大训练组会造成\(J_{train}(\Theta)\)上升,\(J_{CV}(\Theta)\)会持续下降。\(J_{train}(\Theta) < J_{CV}(\Theta)\)然后误差依旧明显存在。
为高方差模型提供更多的训练数据可以解决问题。

神经网络:
小型神经网络(小特征)易于产生高偏差问题,同时计算成本很低。
大型神经网络(多特征)易于产生高方差问题,同时计算成本很高。

机器学习系统设计

偏态类别(skewed classes)的错误分析

如果使用数据组不能涵盖所有/整体的数据组别,偏态类别的情况会发生。也就是说有一个类别的数据远远多于其他类别的数据。

精确率和召回率(Precision and Recall)

预测:1, 实际:1 —— 真正(True Positive)
预测:0, 实际:0 —— 真负(True Negative)
预测:0, 实际:1 —— 假负(False Negative)
预测:1, 实际:0 —— 假正(False Positive)

精确率:
\[\dfrac{\text{True Positives}}{\text{Total number of predicted positives}}
= \dfrac{\text{True Positives}}{\text{True Positives}+\text{False positives}}\]
召回率:
\[\dfrac{\text{True Positives}}{\text{Total number of actual positives}}= \dfrac{\text{True Positives}}{\text{True Positives}+\text{False negatives}}\]
准确率:
\[\frac {true positive + true negative} {total population}\]

精确率和召回率的权衡

提高判断的阀限可以提高精确率,但是会降低召回率。模型的预测确信度提高。
降低判断的阀限可以提高召回率,但是会降低精确率。模型的预测安全性提高。

F数值(也叫F1数值)是用于计算一个模型的精确率和召回率的指标。
\[\text{F Score} = 2\dfrac{PR}{P + R}\]
在验证组上算出F数值可以避免在测试组上出现精确率和召回率的问题。

一个模型需要采取足够的特征才能有足够的信息去判断。常用的方法:如果给予输入X,一个人类专家能否自信地判断出y?
大数据的原则:如果模型的偏差很低,给予越多的数据就越能减少过拟合的问题。

第七周

支持向量机(Support Vector Machine)(SVM)

先设置z:
\[\begin{gather}
z = \theta^Tx \cr
\text{cost0} (z) = \max(0, k(1+z)) \cr
\text{cost1} (z) = \max(0, k(1-z))
\end{gather}
\]
代价函数:
\[J(\theta) = C\sum_{i=1}^m y^{(i)} \ \text{cost1}(\theta^Tx^{(i)}) + (1 - y^{(i)}) \ \text{cost0}(\theta^Tx^{(i)}) + \dfrac{1}{2}\sum_{j=1}^n \Theta^2_j\]
其中\(C = \frac{1}{\lambda}\)。这条公式已经是优化的版本。
假设函数:
\[h_\theta(x) =\begin{cases} 1 & \text{if} \ \Theta^Tx \geq 0 \ 0 & \text{otherwise}\end{cases}\]

大边界分类

如果C的数值极大,代价函数的公式则是:
\[\begin{align*}
J(\theta) = C \cdot 0 + \dfrac{1}{2}\sum_{j=1}^n \Theta^2_j \newline
= \dfrac{1}{2}\sum_{j=1}^n \Theta^2_j
\end{align*}\]

核函数(Kernels)

使用该方法的SVM能制造非线性分类方法。

先使用多个坐标\(l^{(1)},\ l^{(2)},\ l^{(3)}\),计算出新特征和特征之间的邻近值。
\[f_i = similarity(x, l^{(i)}) = \exp(-\dfrac{||x - l^{(i)}||^2}{2\sigma^2})\]
使用高斯核函数,相似度公式则改成:
\[f_i = similarity(x, l^{(i)}) = \exp(-\dfrac{\sum^n_{j=1}(x_j-l_j^{(i)})^2}{2\sigma^2})\]
假设函数就变成:
\[
\begin{align*}l^{(1)} \rightarrow f_1 \newline l^{(2)} \rightarrow f_2 \newline l^{(3)} \rightarrow f_3 \newline\dots \newline h_\Theta(x) = \Theta_1f_1 + \Theta_2f_2 + \Theta_3f_3 + \dots\end{align*}
\]

核函数(Kernels II)

可以把坐标设置在训练数据上,然后每一个训练数据都有一个坐标。\(f_1 = similarity(x,l^{(1)}) ...\)
\[
x^{(i)} \rightarrow \begin{bmatrix}f_1^{(i)} = similarity(x^{(i)}, l^{(1)}) \newline f_2^{(i)} = similarity(x^{(i)}, l^{(2)}) \newline\vdots \newline f_m^{(i)} = similarity(x^{(i)}, l^{(m)}) \newline\end{bmatrix}
\]
用\(f^{(i)}\)来取代\(x^{(i)}\):
\[\min_{\Theta} C \sum_{i=1}^m y^{(i)}\text{cost1}(\Theta^Tf^{(i)}) + (1 - y^{(i)})\text{cost0}(\theta^Tf^{(i)}) + \dfrac{1}{2}\sum_{j=1}^n \Theta^2_j\]

如果C太大,会造成高方差,低偏差。
如果C太小,会造成高偏差,低方差。

要多元分类的话,参考逻辑回归中的方法(one vs all)。

如果n的数值很大(相对于m),选择逻辑回归或者不包含核函数的SVM。
如果n的数值很小(m的数值中等),选择包含高斯核函数的SVM。
如果n的数值很小(m的数值很大),人为增加特征的数量,用逻辑回归或者不包含核函数的SVM。

]]>
<![CDATA[Andrew Ng机器学习课程笔记(第一至第三周)]]>https://codeproducers.com/20171228/5e0456874b43c3452b5136d8Thu, 28 Dec 2017 07:28:05 GMTAndrew Ng机器学习课程笔记(第一至第三周)

本笔记基于Coursera上吴恩达(Andrew Ng)的机器学习课程内容集结而成。(课程链接)
全课程分为十一周,每周都有1-3个测验。第二周至第九周附有编程作业,一共8份。

第一周:

机器学习导论

机器学习定义:一个电脑程序从任务T学习,获得经验E来提升任务T的表现,整个过程用法方法P来衡量。("A computer program is said to learn from experience E with respect to some class of tasks T and performance measure P, if its performance at tasks in T, as measured by P, improves with experience E.")
监督学习(Supervised Learning)大致分为两种:回归和分类。机器学习需要已标签的数据用作模型训练。(已标签的数据是指预先知道输入和所对应输出的数据组。)
非监督学习(Unsupervised Learning)则并不知道训练数据的输出,而就是说,并不预先知道数据的结果。模型要透过数据来寻找数据之间的结构。

单一变量(单元变量)的线性回归(Linear Regression)

模型:输入单一一个变量x,给出相应的单一结果y。
假设函数(Hypothesis Function):
\[ \hat y = h_\theta(x) = \theta_0 + \theta_{1}x \]
该方程将会是一条直线(\( y = a + bx\))。

代价函数(Cost function)

假设函数的准确性将用代价函数来衡量。
\[ J(\theta_0,\theta_1) = \frac{1}{2m} \sum_{i=1}^m (\hat{y_i} - y_i)^2 = \frac{1}{2m}\sum_{i=1}^m (h_{\theta}(x_i)-y_i)^2 \]
在\(\frac{1}{2}\bar{x}\)中,\(\bar{x}\)是指预测数值和实际数值的平均误差,也就是\((h_{\theta}(x_i)-y)^2\)的平均值。而\(\frac{1}{2m}\)是为了方便之后的微积分计算,因为\(\frac{1}{2}\)会在微积分计算中被约除。

梯度下降(Gradient Descent)

梯度下降法是用代价函数的梯度(slope),逐步修改参数。这样我们可以寻找出在最小误差/最大误差(global/local maximum/minimum**)下,参数的数值。
**(注:梯度下降法并不保证每次都能寻找到全域最大\最小值)
梯度下降的公式
\[\theta_j := \theta_j - \alpha\frac{\partial}{\partial\theta_j} J(\theta_0, \theta_1) \]
在单一变量的线性回归的梯度下降法则为:
\[\begin{gather}
\theta_0 := \theta_0 - \alpha \frac{1}{m} \sum\limits_{i=1}^{m}(h_\theta(x_{i}) - y_{i}) \cr
\theta_1 := \theta_1 - \alpha \frac{1}{m} \sum\limits_{i=1}^{m}((h_\theta(x_{i}) - y_{i}) x_{i})
\end{gather}\]
其中\(\alpha\)的数值是自行设定的,详细在后面讲解。梯度下降法将不停循环直至收敛。

第二周:

多变量的线性回归

在这个情况下,假设函数为:
\[h_\theta (x) = \theta_0 + \theta_1 x_1 + \theta_2 x_2 + \theta_3 x_3 + \cdots + \theta_n x_n\]
用矩阵表示:
\[h_\theta(x) =\begin{bmatrix}\theta_0 \hspace{2em} \theta_1 \hspace{2em} ... \hspace{2em} \theta_n\end{bmatrix}\begin{bmatrix}x_0 \newline x_1 \newline \vdots \newline x_n\end{bmatrix}= \theta^T x\]
注:课程中\(x_0\)的数值假设为1.
为了加快计算,推荐使用数组\(X\)来代替单一输入\(x\)的,这样就不用写循环来计算每一个单一\(x\)了。
假设\(X\)的由3个\(x\)组成,每一行代表一个\(x\),每个\(x\)有1个变量。\(\theta\)依旧使用矩阵来表示。
\[X = \begin{bmatrix}x^{(1)}_0 & x^{(1)}_1 \newline x^{(2)}_0 & x^{(2)}_1 \newline x^{(3)}_0 & x^{(3)}_1 \end{bmatrix},\theta = \begin{bmatrix}\theta_0 \newline \theta_1 \newline\end{bmatrix}\]

新方法的计算公式为:
\[h_{\theta}(X)=X\theta\]

代价函数

新公式为:
\[J(\theta) = \frac{1}{2m}\sum_{i=1}^m (h_{\theta}x^{(i)}-y^{(i)})^2\]
使用矩阵计算的版本为:
\[J(\theta) = \frac {1}{2m} (X\theta - \vec{y})^{T} (X\theta - \vec{y})\]
其中,\(\vec{y}\)是所有y数值的向量。

多变量的梯度下降

梯度下降并没有改变,还是需要为每个参数做运算。
\[\begin{gather}
\theta_0 := \theta_0 - \alpha \frac{1}{m} \sum\limits_{i=1}^{m}(h_\theta(x^{(i)}) - y^{(i)}) \cdot x_{0}^{(i)}) \cr
\theta_1 := \theta_1 - \alpha \frac{1}{m} \sum\limits_{i=1}^{m}((h_\theta(x^{(i)}) - y^{(i)}) \cdot x_{1}^{(i)}) \cr
\theta_2 := \theta_2 - \alpha \frac{1}{m} \sum\limits_{i=1}^{m}((h_\theta(x^{(i)}) - y^{(i)}) \cdot x_{2}^{(i)}) \cr
...
\end{gather}\]
也就是说:
\[ \theta_j := \theta_j - \alpha \frac{1}{m} \sum\limits_{i=1}^{m} (h_\theta(x^{(i)}) - y^{(i)}) \cdot x_j^{(i)} \hskip1em \text{for j := 0..n}\]
使用矩阵计算:
\[\theta := \theta - \frac{\alpha}{m} X^{T} (X\theta - \vec{y})\]

特征标准化(Feature Normalization)

如果数据特征的数值都在差不多的范围内的话,梯度下降法的速度将会加快。这是因为\(\theta\)会在小数值范围的时候下降得更快,在大数值范围的时候速度则变慢。所以当特征数值都不均匀的时候,计算效率会明显下降。
理想情况下,特征数值的范围应为:
\[-1 \le x_{(i)} \le 1\]
或者:
\[-0.5 \le x{(i)} \le 0.5\]
有两种方法可以达到这个效果,分别是:特征缩放(feature scaling)和均值归一(mean normalization)。两种方法合并的话,公式为:
\[x_i := \frac{x_i - \mu_i}{s_i}\]
\(\mu_i\)是所有特征(i)的平均数值。\(s_i\)是特征数值的范围(最大值-最小值),或者使用标准偏差(standard deviation)。

有关梯度下降法的提示

画图。\(J(\theta)\)为y轴,循环的次数为x轴,如果\(J(\theta)\)上升,说明\(\alpha\)的数值过大。
自动收敛测试。如果\(J(\theta)\)的下降值少于E,停止循环。E的数值应小于\(10^{-3}\)。

正规方程(Normal Equation)

使用正规方程可以直接计算出\(\theta\)的最优解,这方法并不需要循环。
公式为:
\[\theta = (X^{T} X)^{-1} X^{T} y\]

正规方程和梯度下降的对比

梯度下降 正规方程
需要选择\(\alpha\) 不需要选择\(\alpha\)
需要许多循环 不需要循环
\(O(kn^2)\) \(O(n^3)\),需要计算\((X^{T}X)^{-1})\)
当\(n\)很大时,运算没有问题 当\(n\)很大时,运算变慢

正规方程中,不可逆矩阵(noninvertible)的计算方法
\(X^{T}X\)可能得出一个不可逆的矩阵。通常是因为特征之间线性相关(linear dependent),有些特征是多余的(Redundant)。或者因为矩阵中有过多的特征(e.g. \(m \le n\),特征的数量多于数据的的数量 ),这种情况下需要删除一些特征或者正规化数据(Regularization)。

使用'pinv'比使用'inv'更好。

第三周

逻辑回归

逻辑回归使用来处理分类问题。

二元分类(Binary Classification)

输出的结果y只会出现0或1。\(y \in {0, 1}\)通常0代表否,1代表是。
假设函数的范围:
\[0 \le h_{\theta}(x) \le 1 \]
假设函数是:
\[\begin{gather}
h_{\theta}(x) = g(\theta^{T}x) \cr
z = \theta^{T}x \cr
g(z) = \frac {1} {1 + e^{-z}}
\end{gather}\]
假设函数实际上给的是输出是1的几率。如果\(h_{\theta}(x) = 0.7 \),即是说有70%的几率输出是1。
\[h_{\theta}(x) = P(y = 1|x; \theta) = 1 - P(y = 0|x; \theta)\]
\[P(y = 0|x; \theta) + P(y = 1|x; \theta) = 1 \]

决策边界

因为最终的结果只能是0或1,所以我们需要将\(h_{\theta}(x)\)的结果转换成0或者1。
\[\begin{gather}h_{\theta}(x) \ge 0.5 \rightarrow y = 1 \cr h_{\theta}(x) < 0.5 \rightarrow y = 0 \end{gather}\]
所以,可以这样表示:
\[\begin{gather}\theta^{T}x \ge 0 \Rightarrow y = 1 \cr \theta^{T}x < 0 \Rightarrow y = 0 \end{gather}\]

代价函数

逻辑回归的代价函数是:
\[\begin{gather} J(\theta) = \dfrac{1}{m} \sum_{i=1}^m Cost(h_\theta(x^{(i)}),y^{(i)}) \cr Cost(h_\theta(x),y) = -\log(h_\theta(x)) \hskip1em \text{if y = 1} \cr Cost(h_\theta(x),y) = -\log(1-h_\theta(x)) \hskip1em \text{if y = 0}
\end{gather}\]
将以上两条合并起来的话,得出的结果是:
\[Cost(h_{\theta}(x), y) = - y \space log(h_{\theta(x)}) - (1 - y)log(1 - h_{\theta}(x)) \]
使用矩阵运算的话:
\[\begin{gather}h = g(X\theta) \cr
J(\theta) = \frac{1}{m} \cdot (-y^{T}log(h) - (1 - y)^{T}log(1-h))
\end{gather} \]

梯度下降

梯度下降的公式为:
\[\theta_j := \theta_j - \alpha \dfrac{\partial}{\partial \theta_j}J(\theta)\]
把偏积分的公式扩展:
\[\theta_j := \theta_j - \frac{\alpha}{m} \sum_{i=1}^m (h_\theta(x^{(i)}) - y^{(i)}) x_j^{(i)}\]
使用矩阵运算:
\[\theta := \theta - \frac{\alpha}{m} X^{T} (g(X \theta ) - \vec{y})\]

多元分类

面对多元分类,需要做的只是为每个类别逐个做一次单独分类。
\[\begin{gather} y \in \lbrace0, 1 ... n\rbrace \cr h_\theta^{(0)}(x) = P(y = 0 | x ; \theta) \cr h_\theta^{(1)}(x) = P(y = 1 | x ; \theta) \cr \cdots \cr h_\theta^{(n)}(x) = P(y = n | x ; \theta) \cr \mathrm{prediction} = \max_i( h_\theta ^{(i)}(x) )\cr\end{gather}\]

正规化(Regularization)

欠拟合则是指假设函数不能很好得在现有数据上给出良好的结果。
过拟合问题是指假设函数能在现有数据上给出良好的结果,但是面对新的数据给出的结果却十分不理想。
能解决过拟合问题的方法有两种:

  1. 减少特征的数量
    a. 人手挑选特征
    b. 使用模型选择算法
  2. 正规化 - 在不减少特征的情况下,减少参数\(\theta_j\)的数值

正规化的线性回归

公式为:
\[min_\theta\ \dfrac{1}{2m}\ \left[ \sum_{i=1}^m (h_\theta(x^{(i)}) - y^{(i)})^2 + \lambda\ \sum_{j=1}^n \theta_j^2 \right]\]
梯度下降改为:
\[\begin{gather}
\theta_0 := \theta_0 - \alpha\ \frac{1}{m}\ \sum_{i=1}^m (h_\theta(x^{(i)}) - y^{(i)})x_0^{(i)} \cr
\theta_j := \theta_j - \alpha\ \left[ \left( \frac{1}{m}\ \sum_{i=1}^m (h_\theta(x^{(i)}) - y^{(i)})x_j^{(i)} \right) + \frac{\lambda}{m}\theta_j \right]
\end{gather}\]
抽出分母:
\[\theta_j := \theta_j(1 - \alpha\frac{\lambda}{m}) - \alpha\frac{1}{m}\sum_{i=1}^m(h_\theta(x^{(i)}) - y^{(i)})x_j^{(i)}\]
正规方程:
\begin{align*}& \theta = \left( X^TX + \lambda \cdot L \right)^{-1} X^Ty \newline& \text{where}\ \ L = \begin{bmatrix} 0 & & & & \newline & 1 & & & \newline & & 1 & & \newline & & & \ddots & \newline & & & & 1 \newline\end{bmatrix}\end{align*}
在原先的正规方程中,当\(m \le n\),\(X^{T}X\)是不可逆的。加上\(\lambda \cdot L\)后,\(X^{T}X + \lambda \cdot L\)是可逆的。

正规化的逻辑回归

公式为:
\[J(\theta) = - \frac{1}{m} \sum_{i=1}^m [ y^{(i)}\ \log (h_\theta (x^{(i)})) + (1 - y^{(i)})\ \log (1 - h_\theta(x^{(i)}))] + \frac{\lambda}{2m}\sum_{j=1}^n \theta_j^2\]
注:第二个和,\(\sum_{j=1}^n \theta_j^2\)并不包括\(\theta_0\),(\(\theta\)从0到n一共有n+1个数值),这个和跳过了\(\theta_0\),从1计算到n。
梯度下降为:
\[\begin{gather}
\theta_0 := \theta_0 - \alpha\ \frac{1}{m}\ \sum_{i=1}^m (h_\theta(x^{(i)}) - y^{(i)})x_0^{(i)} \cr
\theta_j := \theta_j - \alpha\ \left[ \left( \frac{1}{m}\ \sum_{i=1}^m (h_\theta(x^{(i)}) - y^{(i)})x_j^{(i)} \right) + \frac{\lambda}{m}\theta_j \right]
\end{gather}\]

]]>
<![CDATA[为Ghost配置HTTPS]]>https://codeproducers.com/20171122/5e0456874b43c3452b5136d7Wed, 22 Nov 2017 13:10:50 GMT为Ghost配置HTTPS

更新:
9/12/2017: 取消注释http端口,http链接自动转跳https

自HTTPS诞生以来,是否采用HTTPS还是停留HTTP的争议就没有中断。网站不采用HTTPS的原因主要分为三点:

  1. HTTPS需要CA(Certificate Authority)签发证书,而各CA签发证书需要的费用不菲。
  2. 一旦采用HTTPS,过去使用HTTP的资源可能无法读取(例如图片)。
  3. HTTPS较HTTP慢。

针对以上几点,现时其实均有很好的解决方案。而且经过数年发展,各大网站也都使用HTTPS,使用HTTPS的好处主要为两点:

  1. 更加安全,HTTPS链接不怕被网络服务商劫持,例如加插弹窗广告。(右下角的广告就是):
    为Ghost配置HTTPS
  2. SEO优化,搜索引擎会为采用HTTPS的网址额外加分,例如Google。

考虑到以上几点,作者建议所有Ghost使用者应在尽早改用HTTPS,免去日后HTTP搬迁HTTPS的麻烦。
Ghost配置HTTPS的流程十分简单,只需要:

  1. 安装certbot
  2. 更改nginx配置
  3. 自动更新CA证书
  4. 更改Ghost的HTTP链接(可选)

1. 安装certbot

以上提及使用HTTPS需要向CA申请SSL证书,SSL证书分为3种等级,等级越高,申请时间和费用越高:

  1. 域名验证(DV)
  2. 组织验证(OV)
  3. 扩展验证(EV)

每个网站应根据自身需求申请证书。如同作者一样,是为Ghost博客申请SSL证书的话,最低级的DV验证即可。
(有关SSL内容,有兴趣请查阅wiki)

为了推广HTTPS,消除申请证书的复杂步骤和金钱,Let’s Encrypt组织于2015年诞生。该组织致力于推广互联网的加密连接,为网站提供免费的SSL证书。该组织的赞助商包括Mozilla, CISCO, Chrome, EFF等互联网和网络安全巨头。

作者正式使用该组织派发的SSL证书实现HTTPS。现时不少CA也提供免费SSL证书服务,读者可自行选择。

而EFF作为Let’s Encrypt的主要赞助,也推出certbot进一步方便网站建立HTTPS链接。过程十分简单:

  1. 登入certbot官网
  2. 选择自身使用的软件和操作系统。(作者使用的是nginx和ubuntu 16.04)
    为Ghost配置HTTPS
  3. 选择后将自动转跳到安装指南。

作者将演示一遍nginx和ubuntu 16.04的安装过程,不同配置过程会有所不同。

$ sudo apt-get update
$ sudo apt-get install software-properties-common
$ sudo add-apt-repository ppa:certbot/certbot
$ sudo apt-get update
$ sudo apt-get install python-certbot-nginx 

这几步只需要跟着做就可,不会发生问题。(可能有些机器已经安装software-properties-common,照做无妨)

2. 更改nginx配置

理论上跟着官网流程,接下去只需要输入

sudo certbot --nginx

即可完成SSL配置,但是作者发现ubuntu 16.04会普遍发生这种问题:

Could not open file: /etc/nginx/sites-enabled/default

经调查,更改nginx的配置文件即可解决问题
在nginx根目录中(作者的nginx根目录为/etc/nginx),将nginx.conf的这行

include /etc/nginx/sites-enabled/*;

改为

include /etc/nginx/sites-enabled/*.*;

然后重新输入

sudo certbot --nginx

完成SSL配置,配置过程中certbot会询问邮箱,使用域名,以及是否强制将HTTP链接转为HTTPS。(邮箱之后会有Email认证,同时建议强制将HTTP转为HTTPS)

之后打开/etc/nginx/sites-enabled/目录下的网址conf文件(xxxx.conf),会多了以下几行:

listen 443 ssl; # managed by Certbot
ssl_certificate /etc/letsencrypt/live/www.codeproducers.com/fullchain.pem; # managed by Certbot
ssl_certificate_key /etc/letsencrypt/live/www.codeproducers.com/privkey.pem; # managed by Certbot
include /etc/letsencrypt/options-ssl-nginx.conf; # managed by Certbot
ssl_dhparam /etc/letsencrypt/ssl-dhparams.pem; # managed by Certbot

if ($scheme != "https") {
    return 301 https://$host$request_uri;
} # managed by Certbot

自此HTTPS配置成功。
(作者著:Certbot只是添加了这几行,之前HTTP的接口并没有关闭。如果读者不需要保留HTTPS接口,需要自行注释前面并将代码移位,如图)(由于输入网址后,游览器会默认使用http链接,一旦关闭http接口,很大机会会出现网页不存在的错误,所以不建议关闭http接口,certbot自带的代码会自动将http连接转https。)
将certbot添加的代码略作修改(可选,作者使用了http2,第二行代码是接收ipv6的访问,需要添加‘ipv6only=on default_server’这个后续,否则有些版本的nginx会产生ipv4和ipv6冲突的问题,访问会无限循环。)

listen 443 ssl http2 default_server; 
listen [::]:443 ssl http2 ipv6only=on default_server;

(作者最终版本的nginx配置,根据读者选择,可能和读者的版本有所不同)
为Ghost配置HTTPS

2.1 使用HTTP2 (可选)

另外,使用HTTPS本身较HTTP慢,建议使用http2取代(原因参考这篇文章)。方法十分简单,只需要添加http2字眼即可:

listen 443 ssl http2; # managed by Certbot
listen [::]:443 ssl http2;

3. 自动更新CA证书

现时网址已经采用HTTPS链接了,但是事情尚未解决,由Let’s Encrypt所发出的证书仅有90日有效期。所以,还需要编写脚本来自动更新证书。
首先可运行这指令测试certbot能否更新证书:

sudo certbot renew --dry-run

出现Congratulations等则证明测试成功,接下来需要编写脚本。两种方法:

  1. 直接修改crontab文件
  2. 建立script然后修改crontab文件

3.1 直接修改crontab文件:

输入command

sudo crontab -e

然后修改crontab文件,加入以下几行:

0 0 * * * certbot renew --force-renewal
5 0 * * * service nginx restart

0 0 * * * 是crontab的时间指令,五位数分别代表分,时,日,月,星期,* 代表每次。
(图片来自 http://www.cnblogs.com/xiaohaillong/p/6102697.html )
为Ghost配置HTTPS
0 0 * * * 就是每日0点0分执行指令,5 0 * * * 则是每日0点5分执行指令。(预留5分钟时间更新证书)

保存修改就大功告成。

3.2 建立script然后修改crontab文件

在/root目录建立新文件,命名为letsencrypt.sh,文件内容为:

#!/bin/bash
systemctl reload nginx

# 如果有其他进程也使用该证书:
# systemctl restart mosquitto

保存文件,然后更改该文件的使用权限:

chmod +x /root/letsencrypt.sh

输入指令:

sudo crontab -e

修改crontab文件,加入一行:

20 3 * * * certbot renew --noninteractive --renew-hook /root/letsencrypt.sh

(每日3点20分更新证书,然后执行script文件重启nginx)

保存修改就大功告成。

4. 更改Ghost的HTTP链接(可选)

(可选)自此GHOST已换上HTTPS,但不排除在ghost.conf上使用的还是HTTP的链接。
更改为(作者著:这方法不会更改ghost的@blog.url,url还是原来的HTTP链接。由于ghost内部link都是根据blog.url做扩展,很大几率会从https连去http的链接,所以不建议使用):
为Ghost配置HTTPS
或者改为(建议):
为Ghost配置HTTPS
之后后请记得更改Design->Navigation的链接。否则所有HTTP链接都不可用,会返回Page Not Found的错误。

参考:

  1. Nginx on Ubuntu 16.04 (xenial)(certbot官方教程)
  2. 如何免费的让网站启用HTTPS
  3. 1125使用的命令行
  4. https是什么?使用https的好处与不足?
  5. How to setup Let's Encrypt for Nginx on Ubuntu 16.04 (including IPv6, HTTP/2 and A+ SLL rating)
  6. Let’s Encrypt
]]>
<![CDATA[Atom 编辑器插件推荐]]>https://codeproducers.com/20171116/5e0456874b43c3452b5136d6Thu, 16 Nov 2017 08:35:34 GMTAtom 编辑器插件推荐

Atom作为一款开源编辑器,被视为Emacs和Sublime的强有力竞争者。现时经过多年发展,Atom的社区不断壮大,最为明显的就是Atom拥有数以千记的插件。在这些插件当中,作者将推荐自己现时正在使用的插件,希望能帮助读者配置自己适用的Atom运行环境。

注意:

  1. 作者现时主要使用Python和C语言,其他语言使用情况较少,插件使用范围有一定局限性。
  2. 推荐插件标准纯属为作者喜好,难免出现偏颇情况,请读者自行决定是否采用该插件。
  3. 由于Atom更新,插件没有更新等原因,推荐插件可能会出现无法使用等情况,请留意。

实用类 :

  1. atom-beautify
    简介: 数十种语言的代码规整插件。
    优点: Atom居家旅行必备,代码格式规整一click到位。
    缺点: 原生软件仅内置部分语言规格,其他语言需要额外安装。(官网有教程)
    注意: 该软件可能会引发大括号是否新开一行,space还是tab,什么时候空一格等战争风险,请酌量使用。

  2. atom-ide-ui
    简介: 界面功能,代码分析,后台功能等N合一界面插件。
    优点: 十分强大的ide-ui整合,由facebook人马维护,值得信赖。
    缺点: 现时能支持的语言不多。部分语言有第三方支援,需要额外安装。
    注意: 不兼容linter套件,只能二选一。安装完毕会提示停止使用linter,两者一起运行会引发冲突。

  3. auto-encoding
    简介: 读取日文文档,代码必备。Atom虽然内置Shift-JIS转码,但有些日语文字仍然不能正确地显示。安装这个插件能解决问题。
    优点: 存在即是优点。
    缺点: 不需要的话即是缺点。

  4. autocomplete-python
    简介: python语法,语句自动补全插件。
    优点: python码农必备,懒人养成插件。
    缺点: 自动识别有时过于强大,容易误触。(随便打个字母都可以出一堆补全选项)

  5. docblockr
    简介: 代码文件的文档,补充说明,注释补全插件。
    优点: 码农必备,懒人养成插件II。(然而懒人根本不会写注释)
    缺点: 支援的语言不多,暂时只支援主流语言。

  6. Git Plus
    简介: 单靠Atom命令实现git的操作插件。
    优点: 纯CLI操作
    缺点: 纯CLI操作,不习惯CLI的用户可能无所适从。(这里不介入CLI和GUI的战争,人各有好,各取所需)

  7. Highlight Selected
    简介: Highlight插件,一旦出现和选中的字眼一样的字,会一起highlight
    优点: Atom居家旅行必备。(其实是原生Atom太朴素,什么都没有)
    缺点: 不需要的话即是缺点。
    注意: 建议配合一起安装Minimap和minimap-highlight-selected,效果更佳。

  8. Language Haskell
    简介: Haskell的highlight插件
    优点: 写Haskell的话就只有这款靠谱。插件内容可以自定义,官网有教程。
    缺点: 不需要的话即是缺点。

  9. language-scheme
    简介: Scheme的highlight插件
    优点: Scheme的highlight好像就只有这款插件。(16/11/2017:数了一数,Atom上有关Scheme的插件好像只有个位数。)
    缺点: Scheme插件都没几个,不能要求这么多。(Make Scheme Great Again)

  10. linter
    简介: 界面功能,代码分析,后台功能等N合一界面插件。
    优点: Atom很早就存在的界面插件,使用人数众多,第三方配套插件齐全。
    缺点: linter的第三方配套插件并没有那么好用,设置比较麻烦。
    注意: 安装后会询问用户是否一同安装Busy Signal,Intentions,Linter Ui Default三款插件(建议安装)。不兼容atom-ide-ui套件,只能二选一.安装完毕会提示停止使用atom-ide-ui,两者一起运行会引发冲突。

  11. Markdown Preview Enhanced
    简介: Markdown预览插件。
    优点: Markdown支援齐全。对比几款Atom上的Markdown插件,作者用这款的原因更多是看得舒服。(比原生Atom那破预览不知道高到哪里去了)
    缺点: 不需要的话即是缺点。
    注意: 不兼容原生Atom的Markdown Preview插件,安装后需要Disable原生Atom插件。

  12. Markdown-Writer
    简介: Markdown语法补全。
    优点: 码农必备,懒人养成插件III。
    缺点: 不需要的话即是缺点。作者装了之后其实没怎么用。

  13. Minimap
    简介: Minimap显示。
    优点: 存在即是优点。可设置,自定义的地方多。
    缺点: 不需要的话即是缺点。
    注意: 建议一起安装Highlight Selected和minimap-highlight-selected,效果更佳。

  14. minimap-highlight-selected
    简介: Minimap显示支援highlight。
    优点: 存在即是优点。
    缺点: 不需要的话即是缺点。
    注意: 安装前需要安装Highlight Selected和Minimap才能生效。

  15. pigments
    简介: 令颜色代码显示出相对应的颜色。
    优点: 写CSS必备,其他语言写RGB代码一样支援。自定义的选择多。
    缺点: 默认配置有时会将不关乎颜色的代码一样着色。例如python的注释,用#开头的注释容易出现这种情况。

  16. PlatformIO IDE Terminal
    简介: 顾名思义,将terminal界面整合到Atom。
    优点: 一键安装,容易使用。自定义的选择多。
    缺点: 不需要的话即是缺点。
    注意: 请勿使用terminal-plus这款插件,作者已经数年没有更新,并全力投入到PlatformIO IDE Terminal的研发上。

  17. Project Manager
    简介: 项目管理神器。
    优点: 多个项目轻松切换,容易使用管理。自定义选择多。
    缺点: 切换项目之后需要等Atom缓慢加载。(Atom的锅)

  18. Script
    简介: 代码运行插件,无需terminal加载。
    优点: 直接界面运行,不要另外开terminal。支援的语言众多。
    缺点: 原装插件内置的支援语言不多,需要安装第三方插件。很多语言只能在MacOS或Linux运行环境下运行,详情参考官网。

  19. Todo Show
    简介: 注释中的TODO,FIXME,BUG等字眼管理插件。
    优点: 一键寻找前人/自己N年前留下的锅注释,哪里地方需要更改一目了然。
    缺点: 有了这款软件人类还是会留下TODO等字眼不fix。一些library库中的字眼也会一同显示,需要额外设置。

装饰类(没啥用类):

  1. Activate Power Mode
    简介: 纯粹好玩没用的软件,打字combo插件。
    优点: 哇,coding好有成就感哦,看我python1000连击(大误)。
    缺点: Atom打字卡机lag机的万恶之源。建议鶸电脑/笔电关闭震动和特效,保留combo数即可。

  2. Atom Clock
    简介: 界面上的时间显示插件。
    优点: 可自定义。其实看着还挺萌。
    缺点: 工作时程序员根本就不会留意时间。(文明玩家:duang,下一回合。程序员:duang,下一个bug。)

  3. File Icons
    简介: 图标补充插件。
    优点: 图标应有尽有,基本上所有文件后续都有对应的图标。这类别最有用的插件。
    缺点: 然而还是没什么大作用,纯粹为了好看。

参考:

  1. Atom编辑器
  2. Atomで楽しくC,C++開発をする個人的設定まとめ
  3. Atom で Markdown を 快適に書くためのTips
  4. Atom Packages Most Stars
]]>
<![CDATA[使用ghost建立一个免费博客网站]]>https://codeproducers.com/20171108/5e0456874b43c3452b5136d5Wed, 08 Nov 2017 09:06:03 GMT使用ghost建立一个免费博客网站

现时建立一个博客的方法五花八门,大致可分为三类。

  1. 使用第三方博客平台,直接申请一个博客网站,例如Blogger, Medium
  2. 使用第三方博客框架,建立一个静态博客网站,可以将Github内容转化为博客网站,例如jekyll, GitHub Pages
  3. 使用第三方网站框架,建立一个博客网站,可以根据自己需要设置网站功能,例如WordPress, ghost

本文章采用的是第三种方法,使用ghost框架来搭建一个免费博客网站。

为什么使用ghost

在以上介绍的三类方法中,以第一种最为方便,注册账号即可建立自己的博客网站,但可修改的空间不多。第三种最为自由度最高,但是建立网站的步骤比较繁琐。第二种为两方面的折中选择。读者可按照自身情况建立自己的博客。如读者自身有一定编程基础,建议使用第二种或第三种方法。

在各种第三方网站框架中,ghost相比其他框架,优势在于自身的小体积和易于上手。相比WordPress那臃肿的框架,ghost的目标十分明确,就是建立一个轻量级,易于上手又不乏扩展性的博客平台。

建立一个博客网站的流程

使用ghost建立一个网站大致分为以下几个步骤,部分流程的先后次序并非绝对,读者可按自身需求调整。(方法一)

  1. 为博客申请一个(免费)伺服器(可选)
  2. 在伺服器上安装必要前置软件
  3. 在伺服器上搭建ghost平台
  4. 为博客购买和设置域名(可选)

或者可以选用Bitnami软件包,流程为:(方法二)

  1. 为博客申请一个(免费)伺服器(可选)
  2. 安装Bitnami软件包
  3. 为博客购买和设置域名(可选)

虽然Bitnami简化了安装流程,实现了一键安装ghost,但是作者在使用Bitnami之后,出现ghost一旦关闭就不能重新开启的bug(Issue 501),所以不推荐使用Bitnami。

1. 申请一个(免费)伺服器(可选)

博客平台可以建立在自己的电脑或者搭建在第三方伺服器上,两者各有优势。

搭建在自家电脑的优势:

  • 不需要远程管理,伺服器的情况自己可以在地监控。

搭建在自家电脑的缺点:

  • 需要24小时开机运行,一旦关机,博客就不能被访问。

搭建在第三方平台的优势:

  • 第三方平台能提供强大的技术支援和功能扩展。

搭建在第三方平台的缺点:

  • 通常需要一定费用来维持伺服器运作
  • 需要远程管理

本文选择使用亚马逊网上服务(AWS)的EC2来搭建一个免费的伺服器(中国读者请使用中国亚马逊云服务)。现时AWS推出不少免费的服务,其中包括了EC2,免费期为期一年。读者可以按照自己喜好使用其他伺服器平台,例如阿里云等。

使用AWS服务需要:

  1. 常用邮箱
  2. 信用卡
  3. 电话号码

1.1 注册一个AWS服务

登入AWS,注册一个免费账号。
使用ghost建立一个免费博客网站
之后填写个人账号名称,邮箱,密码。进入下一页。
使用ghost建立一个免费博客网站
选择个人账号,填写个人的联络资讯。
使用ghost建立一个免费博客网站
填写个人信用卡资料,AWS会向信用卡账号收取1美元做验证,之后会归还。
使用ghost建立一个免费博客网站
然后是电话验证,AWS会打电话到填写的电话号码,这时候需要输入AWS网上提供的验证码。(图片来自apachefriends)
使用ghost建立一个免费博客网站
接着选择Basic服务。(图片来自apachefriends)
使用ghost建立一个免费博客网站
最后确认所有资讯无误,注册AWS账号成功。

1.2 申请一个EC2伺服器

登入AWS Console,选择EC2。
使用ghost建立一个免费博客网站
选择'Launch Instance',建立一个新伺服器。
使用ghost建立一个免费博客网站

现时可以建立两种伺服器:

  1. 如果不选用Bitnami软件包,建立一个ubuntu伺服器。(方法一)
  2. 如果选用Bitnami软件包,建立一个自带ghost的伺服器。(方法二)(不建议)
方法一:建立ubuntu伺服器

选择ubuntu 16.04
使用ghost建立一个免费博客网站
注意:请谨记选择t2.micro,唯有t2.micro才能享有AWS一年免费服务,即使选择t2.nano也会收钱。
使用ghost建立一个免费博客网站
然后不要选'Review and Launch',先设置其他。'Configure Instance'可以不用理会,'Add Storage'则是选择伺服器容量,不建议太小,也不建议太大。建议8GB-10GB即可。
注意:免费服务只提供30GB容量,请不要超出。
使用ghost建立一个免费博客网站
'Configure Security Group'方面,添加HTTP和HTTPS接口,否则伺服器不能透过输入网址的方式访问。
使用ghost建立一个免费博客网站
确认无误后,点右下角的launch便可建立server。
使用ghost建立一个免费博客网站
最后AWS会询问关于配对钥匙的方案,这是关于SSH的访问使用。原因请参考RSA加密。请选择'Create a new key pair',产生一堆新的钥匙。并为钥匙命名。这时候游览器会自动下载自己的密匙。
注意:请务必保存好密匙文件,因为AWS一旦发出密匙文件,就不会保存。一旦密匙丢失,是没有方法找回来的。而密匙是透过SSH连上伺服器的关键,请妥善保管。
使用ghost建立一个免费博客网站
点击'Launch Instance',ubuntu伺服器就成功建立了。

方法二:使用bitnami,建立自带ghost的伺服器

先点击AWS Marketplace。
使用ghost建立一个免费博客网站
然后在搜寻栏输入ghost,选择'Ghost Certified by Bitnami',
使用ghost建立一个免费博客网站
按continue,
使用ghost建立一个免费博客网站
同方法一,选择t2.micro。其他选项基本上Bitnami已经预先设置好。(Bitnami预设容量为10GB)
使用ghost建立一个免费博客网站
最后AWS同样会询问关于配对钥匙的方案。选择'Create a new key pair',保管好key即可。
使用ghost建立一个免费博客网站
点击'Launch Instance',伺服器就成功建立了。

使用方法二的话,bitnami软件包自带安装apache,mysql,node.js和ghost。一旦伺服器运行会自动运行ghost,直接完成网站建设。

1.3 访问EC2伺服器

建立伺服器后,可透过SSH方法连访问伺服器。支持SSH链接的工具有很多,例如PuTTY或者WinSCP。读者可按照自身需要使用。作者使用PuTTY作为示范。

返回EC2界面,'Running Instances'显示有多少server正在运行。点击'Running Instances'。
使用ghost建立一个免费博客网站
点击server,选择左上'connect',AWS会教如何使用putty登入伺服器。保存Public DNS网址(xxxxxx.xxxxx.compute.amazonaws.com)
使用ghost建立一个免费博客网站

注意:如果使用Bitnami软件包,这里有额外一步-取得应用密码(Application Password),请务必在登入前完成,应用密码仅在第一次显示。
取得应用密码:

  1. 点击'Actions',选'Instance Setting',点'Get System Log'(图片来自@robmadd3n)
    使用ghost建立一个免费博客网站
  2. 在log文件中寻找'Setting Bitnami application password to'(图片来自@robmadd3n)
    使用ghost建立一个免费博客网站
  3. 保存该密码,每次登入都需要输入该密码。

使用PuTTY登入伺服器分为2步:

  1. 将之前下载的密匙(xxx.pem)转换为putty密匙格式(xxx.ppk)
  2. 透过PuTTY,设置伺服器网址和链接密匙,登入伺服器

(步骤可参考AWS文档)

1.3.1 转换密匙

开启PuTTY Key Generator(程序缩写为PuTTYgen.exe)
使用ghost建立一个免费博客网站
点击'Load',选取之前下载的pem文件,PuTTY默认打开.ppk文件,请更改为All Files
使用ghost建立一个免费博客网站
之后PuTTY会跳出提示,提醒要保存为'private key'。跟从提示,点击'Save private key',完成步骤。

1.3.2 PuTTY登入设置

打开PuTTY,输入Public DNS网址在Host Name一栏。
使用ghost建立一个免费博客网站
在左手边的'Category'中,点'Connection',点'SSH',点'Auth'。右手边选'Browse...',打开刚生成的ppk文件。
使用ghost建立一个免费博客网站
返回一开始的PuTTY页面('Session'), 点击Open即可链接上伺服器。(登入前可以保存设置,方便下次登入)
使用ghost建立一个免费博客网站

如果不使用Bitnami(方法一),用户名就是ubuntu;如果使用了Bitnami软件包(方法二),用户名就是bitnami,还需要输入之前保存的应用密码。

Ubuntu登入成功:
使用ghost建立一个免费博客网站

Bitnami用户(方法二)到这一步可以直接调至域名购买和设置,ubuntu用户继续(方法一)。

2. 在伺服器上安装必要前置软件

(可参考ghost官方文档)
注意:作者安装的是production版本,如果需要安装development版本,请阅读ghost官方文档。

建议新用户登入ubuntu伺服器后,先运行一下指令更新系统。

sudo apt-get update && sudo apt-get dist-upgrade

在安装ghost前,需要安装的软件为node.js,mysql,nginx

2.1 安装node.js

安装node.js十分简单,直接输入命令

curl -sL https://deb.nodesource.com/setup_6.x | sudo -E bash 

再输入命令,直接安装node.js

sudo apt-get install -y nodejs

注意:ghost官方不建议使用nvm来安装,更新和维护node.js,建议删除nvm以免造成影响。详情参考官方文档。

2.2 安装mysql

安装mysql只需要输入一项命令即可

sudo apt-get install mysql-server

过程中,系统会询问用户,要求设置一个密码,请谨记该密码。

2.3 安装nginx

如同安装mysql和node.js,输入命令

sudo apt-get install nginx

以免防火墙阻挡外部访问链接,输入命令

sudo ufw allow 'Nginx Full'

自此所有前置已经安装完毕,下一步可以安装ghost

3. 在伺服器上安装ghost

先安装ghost-cli,这会方便输入ghost相关指令

sudo npm i -g ghost-cli

建议新开一个文件夹以存放ghost和相关插件,输入命令,将会建立一个名加ghost的文件夹在/var/www的目录下。(如果/var/www不存在,则会一同建立。)

sudo mkdir -p /var/www/ghost

注意:ghost不能安装在server根目录(即是/root),安装会出现bug,不能安装等问题。

移至新文件夹

cd /var/www/ghost

输入ghost-cli命令,安装ghost

ghost install

安装过程中会要求输入一下资料(ghost版本的差异可能导致问题顺序不同,新版本会要求回答更多问题)
使用ghost建立一个免费博客网站(图片来自ghost官方)
(过程可参考ghost-install文档)
首先是博客域名'Enter your blog Url:',如果暂时没有想好或申请,可以只打ip地址暂代。
(ghost可能会询问'Do you wish to set up ssl?', 输入'no')
其次是MySQL设置'Enter your MySQL hostname [localhost]:',如果按照以上流程,MySQL安装在本地伺服器,直接按'Enter'键即可。
接着是'Enter your MySQL username:',按照以上流程,输入'root'。如果之前有为数据开新开用户,则使用该用户名。
跟着是'Enter your MySQL password:',输入数据库密码。
再者是'Ghost database name:',为ghost的database起名。
之后是'Do you wish to set up a ghost MySQL user?',为ghost建立一个新用户。
然后是'Do you wish to set up nginx?',这里可以输入yes,之后也可以更改。
'Do you wish to set up systemd?'这个可以输入yes,之后也可以更改。
'Do you want to start Ghost?'输入yes,运行ghost。

以上资料可以在config.development.json更改(位置在ghost安装目录下)
(参考ghost-config文档)

这时候输入伺服器ip:ghost端口就能正常访问博客。

如出现任何问题,参考ghost-troubleshooting文档。

4. 为博客购买和设置域名(可选)

域名可以在网上各大域名商购买,例如godaddy,namecheap等。作者的域名是在namesilo购买。在不同平台购买域名都各有优势,读者可按照自身需求购买。各大域名平台都长期有优惠卷提供,能便宜1美金不等,请自行搜索。
注意:在中国内地建立的伺服器购买域名时需要域名备案,外国第三放域名平台可能不设有域名备案服务,请在购买域名前做好搜索,选择合适的域名商。

这里以namesilo为例,注册账号和购买域名(过程略),登入后选择'domain manager'。
使用ghost建立一个免费博客网站
点选右边第三个按钮(Manage DNS for this domain)
使用ghost建立一个免费博客网站
namesilo会为所有网址预设广告link等,这里全部删除。
然后选取'A'
使用ghost建立一个免费博客网站
输入AWS EC2的伺服器ip(仅限ipv4地址)
使用ghost建立一个免费博客网站
建议一同开启www网址前置
使用ghost建立一个免费博客网站

注意:AWS EC2中每次伺服器重启都会委派新的ip地址,如果想ip地址固定,可以选用EC2的elastic IP,令ip地址固定,不用更新DNS。(教程略)

如果之前ghost url设置上使用的是ip地址,请更改ghost-config文件和nginx配置。(config.development.json位置在ghost安装目录下,nginx配置ghost.conf在/etc/nginx/sites-available/)
(nginx配置可参考ghost中国)

DNS更新需要一定时间,如果半小时至一小时后,输入域名仍未能访问ghost,则很有可能是设置上有问题,建议重新检视DNS,nginx和ghost的设定。

如有任何疑惑和问题,欢迎在讨论区讨论。

参考

  1. Register to Use Amazon Web Services (AWS) and Create an Amazon Account
  2. Host your Application in the Amazon Cloud with XAMPP and Bitnami
  3. 基于 Amazon AWS EC2 部署 Ghost 博客
  4. 在亚马逊云服务AWS上使用 EC2 和 RDS 来搭建ghost博客系统支持Markdown编辑
  5. The Ghost Guide
  6. Create a Professional Blog with AWS and Ghost
  7. docs.ghost.org/docs
]]>