一個小知識：

有意思的是，scipy囊括了numpy的命名空間，也就是說所有np.func都可以通過sp.func等價調用。

簡介：

本部分對一個互聯網公司的流量進行擬合處理，學習最基本的機器學習應用。

導入包&路徑設置：

import os
import scipy as sp
import matplotlib.pyplot as plt

data_dir = os.path.join(
    os.path.dirname(os.path.realpath(__file__)), "..", "data")
# __file__ 是用來獲得模塊所在的路徑的，這可能得到的是一個相對路徑，
# os.path.dirname(__file__) ，相對路徑時返回值為空，
# 為了得到絕對路徑，需要 os.path.dirname(os.path.realpath(__file__)）。
# .realpath得到完整路徑加文件名
# .dirname會去掉腳本名，只保留路徑
# print(os.path.realpath(__file__))
# print(os.path.dirname(os.path.realpath(__file__)))

讀入&清洗數據：

這裏沒有采用原書的讀取方式，似乎原作者不知道這樣寫更為簡潔：

# data = sp.genfromtxt(os.path.join(data_dir, "web_traffic.tsv"), delimiter="\t")
# x = data[:,0]
# y = data[:,1]

# 和上面三行等價
x,y = sp.loadtxt(os.path.join(data_dir, "web_traffic.tsv"), delimiter="\t",unpack=True)

print("非法數據：", sp.sum(sp.isnan(y)))    # 統計nan缺失
x = x[~sp.isnan(y)]                        # 布爾索引
y = y[~sp.isnan(y)]                        # 布爾索引

繪圖函數：

這個函數寫的很精妙，有不少使用了python高級技巧的地方，值得學習：

list+zip的運用，list生成器的運用等等

順便一提，sp.poly1d()生成對象有屬性f.order，可以查看自身的階數。

colors = [‘g‘, ‘k‘, ‘b‘, ‘m‘, ‘r‘] #<-----------
linestyles = [‘-‘, ‘-.‘, ‘--‘, ‘:‘, ‘-‘] #<-----------

def plot_models(x, y, models, fname, mx=None, ymax=None, xmin=None):
    ‘‘‘
    繪制原數據散點圖和擬合線圖
    :param x:        橫坐標 
    :param y:        縱坐標
    :param models:   擬合線（list傳入）
    :param fname:    保存圖像名
    :param mx:       擬合線x的list是否給定了 
    :param ymax:     y軸上限
    :param xmin:     x軸下限
    :return:         None
    ‘‘‘
    plt.clf()                                                         # 清空當前坐標上圖像
    plt.scatter(x, y, s=10, alpha=1, marker=‘.‘)
    # c:散點的顏色
    # s：散點的大小
    # alpha:是透明程度
    # plt.title("上個月網絡流量圖")
    plt.xlabel("Time")
    plt.ylabel("Hits/hour")
    plt.xticks(
        [w*7*24 for w in range(10)], ["week %i" % w for w in range(10)])

    if models:                                                       # 是不是繪制擬合線
        if mx is None:
            mx = sp.linspace(0, x[-1], 1000)
        for model, style, color in zip(models, linestyles, colors): #<-----------
            plt.plot(mx, model(mx), linestyle=style, linewidth=2, c=color) #<-----------

        plt.legend(["d=%i" % m.order for m in models], loc="upper left") #<-----------
    plt.autoscale(tight=True)
    plt.ylim(ymin=0)
    if ymax:
        plt.ylim(ymax=ymax)
    if xmin:
        plt.xlim(xmin=xmin)
    plt.grid(True, linestyle=‘-‘, color=‘0.75‘) #<-----------        # 網格線設置，color應該是灰度
    plt.savefig(fname)

全數據繪圖：

# 繪制原始數據散點圖
plot_models(x, y, None, os.path.join( "..", "1400_01_01.png"))

fp1, res, rank, sv, rcond = sp.polyfit(x, y, 1, full=True)
print("擬合參數： %s" % fp1)
print("誤差數值： %s" % res)
f1 = sp.poly1d(fp1)
f2 = sp.poly1d(sp.polyfit(x, y, 2))
f3 = sp.poly1d(sp.polyfit(x, y, 3))
f10 = sp.poly1d(sp.polyfit(x, y, 10))
f100 = sp.poly1d(sp.polyfit(x, y, 100))

plot_models(x, y, [f1], os.path.join("..", "1400_01_02.png"))
plot_models(x, y, [f1, f2], os.path.join("..", "1400_01_03.png"))
plot_models(x, y, [f1, f2, f3, f10, f100], os.path.join("..", "1400_01_04.png"))

轉折點處理：

分離轉折點前後的數據：

‘‘‘
轉折點處理部分
‘‘‘
inflection = 3.5*7*24
xa = x[:inflection]
ya = y[:inflection]
xb = x[inflection:]
yb = y[inflection:]       # 註意切片的寫法，沒有逗號

轉折點前後分批處理：

# 轉折點前一階擬合
fa = sp.poly1d(sp.polyfit(xa, ya, 1))
# 轉折點後一階擬合
fb = sp.poly1d(sp.polyfit(xb, yb, 1))
plot_models(x, y, [fa, fb], os.path.join("..", "1400_01_05.png"))

# 平方差
def error(f, x, y):
    return sp.sum(sp.sum(f(x) - y)**2) #<-----------

print("全數據點誤差統計：")
for f in [f1, f2, f3, f10, f100]:
    print("Error d=%i: %f" % (f.order, error(f, x, y)))

print("轉折點後誤差統計：")
for f in [f1, f2, f3, f10, f100]:
    print("Error d=%i: %f" % (f.order, error(f, xb, yb)))

print("一階拼接擬合誤差統計： %f" % (error(fa, xa, ya) + error(fb, xb, yb)))

‘‘‘
趨勢預測部分
‘‘‘
# 全數據6周預測
plot_models(x, y, [f1, f2, f3, f10, f100], os.path.join("..", "1400_01_06.png"),
            mx=sp.linspace(0 , 6*7*24, 100),
            ymax=10000, xmin=0)
# 全模型轉折點後擬合
fb1 = fb
fb2 = sp.poly1d(sp.polyfit(xb, yb, 2))
fb3 = sp.poly1d(sp.polyfit(xb, yb, 3))
fb10 = sp.poly1d(sp.polyfit(xb ,yb, 10))
fb100 = sp.poly1d(sp.polyfit(xb, yb, 100))
print("全模型轉折點後誤差統計：")
for f in [fb1, fb2, fb3, fb10, fb100]:
    print("Error d=%i： %f" % (f.order, error(f, xb, yb)))

# 轉折點後數據6周預測
plot_models(
    x, y, [fb1, fb2, fb3, fb10, fb100], os.path.join("..", "1400_01_07.png"),
    mx=sp.linspace(0, 6*7*24, 100),
    ymax=10000, xmin=0)

折點後面的數據單獨處理：

sp.random.permutation()這個函數返回打亂的input

import scipy as sp
sp.random.permutation([1,2,3,4,5])
# Out[3]: 
# array([4, 5, 1, 2, 3])

30%用於測試，70%用於擬合，這裏隨機分離數據

‘‘‘
轉折點後數據
部分用於訓練
部分用於測試
‘‘‘
frac = 0.3
split_idx = int(frac*len(xb))                             # 30%的數據量
shuffled = sp.random.permutation(list(range(len(xb))))    # 全xb的index亂序
test = sorted(shuffled[:split_idx])                       # 亂序index提取前30%，後排序
train = sorted(shuffled[split_idx:])                      # 亂序index提取後70%。後排序

擬合對比：

fbt1 = sp.poly1d(sp.polyfit(xb[train], yb[train], 1))
fbt2 = sp.poly1d(sp.polyfit(xb[train], yb[train], 2))
fbt3 = sp.poly1d(sp.polyfit(xb[train], yb[train], 3))
fbt10 = sp.poly1d(sp.polyfit(xb[train], yb[train], 10))
fbt100 = sp.poly1d(sp.polyfit(xb[train], yb[train], 100))

print("測試點誤差：")
for f in [fbt1, fbt2, fbt3, fbt10, fbt100]:
    print("Error d=%i： %f" % (f.order, error(f, xb[test], yb[test])))
# 繪制部分訓練模型擬合圖
plot_models(x, y, [fbt1, fbt2, fbt3, fbt10, fbt100], os.path.join(‘..‘, ‘1400_01_08.png‘),
            mx=sp.linspace(0, 6*7*24),
            ymax=10000, xmin=0)

優化器解方程：

from scipy.optimize import fsolve
# 想要預測訪問量100000的時間
print(fbt2)
print(fbt2-100000)
reached_max = fsolve(fbt2 - 100000, 800) / (7*24)
print("100,000 hits/hour excpeted at week %f" % reached_max)

『Python』MachineLearning機器學習入門_極小的機器學習應用