### ------------------------------------------------------------------- ###
#     利用 Python 練習資料探勘 URL：http://python.jobbole.com/83563/

# 資料匯入和視覺化
import urllib2
url = "http://aima.cs.berkeley.edu/data/iris.csv"
u = urllib2.urlopen(url)
localFiel = open("iris.csv", "w")
localFiel.write(u.read())
localFiel.close()

"""
資料來源： 伯克利大學
資料包含鳶尾花（iris）資料集，包含了三種鳶尾花（山鳶尾、維吉尼亞鳶尾和變色鳶尾）的各 50 個數據樣本的多元資料集
每個樣本有四個特徵，即花萼（sepal）和花瓣（petal）的長度和寬度，以釐米為單位
資料集有 5 列，前四列包含著特徵值，最後一列代表著樣本型別
"""

# csv 檔案很容易被 numpy 庫的 genfromtxt 方法解析
from numpy import genfromtxt,zeros
# 讀取前 4 列
data = genfromtxt("iris.csv", delimiter=",", usecols=(0, 1, 2, 3))
# 讀取第 5 列
target = genfromtxt("iris.csv", delimiter=",", usecols=(4), dtype=str)

print "data: ", data
print "target: ", target
print "data.shape: ", data.shape
print "target.shape: ", target.shape

# 檢視有多少種樣本型別以及它們的名字
print set(target)

# 使用 pylab 庫（matplotlib的介面）的 plotting 方法可以建一個二維散點圖讓我們在兩個維度上分析資料集的兩個特徵值
from pylab import plot, show, close
# 藍色點代表山鳶尾、紅色點代表變色鳶尾、綠色點代表維吉尼亞鳶尾
# 第一和第三維度是花萼的長度和花瓣的長度
plot(data[target == "setosa", 0], data[target == "setosa", 2], "bo")
plot(data[target == "versicolor", 0], data[target == "versicolor", 2], "ro")
plot(data[target == "virginica", 0], data[target == "virginica", 2], "go")
# show()
close()

# 另一種常用的檢視資料的方法是分特性繪製直方圖
# 下面的程式碼可以繪製資料中每一型別的第一個特性（花萼的長度）
from pylab import figure, subplot, hist, xlim, show
xmin = min(data[:, 0])
xmax = max(data[:, 0])
figure()
subplot(411)
hist(data[target == "setosa", 0], color="b", alpha=.7)
xlim(xmin, xmax)
subplot(412)
hist(data[target == "versicolor", 0], color="r", alpha=.7)
xlim(xmin, xmax)
subplot(413)
hist(data[target == "virginica", 0], color="g", alpha=.7)
xlim(xmin, xmax)
subplot(414)
hist(data[:, 0], color="y", alpha=.7)
xlim(xmin, xmax)
# show()
close()

# 分類
#     ------  高斯樸素貝葉斯分類 ------
# 把字串陣列轉型成整型資料
t = zeros(len(target))
t[target == "setosa"] = 1
t[target == "versicolor"] = 2
t[target == "virginica"] = 3

# 模型例項化和訓練分類器
from sklearn.naive_bayes import GaussianNB
classifier = GaussianNB()
classifier.fit(data, t)

# 分類器可以由 predict 方法完成，並且只要輸出一個樣例就可以很簡單的檢測
print classifier.predict(data[0])
print t[0]

# 評估分類器
# 通過從源資料集中隨機抽取樣本把資料分為訓練集和測試集，然後使用訓練集的資料來訓練分類器，並使用測試集來測試分類器
from sklearn import cross_validation
train, test, t_train, t_test = cross_validation.train_test_split(data, t, test_size=0.4, random_state=0)

# 訓練分類器並輸出精確度
classifier.fit(train, t_train)
print classifier.score(test, t_test)

# 另一個估計分類器表現的工具叫做混淆矩陣。在此矩陣中每列代表一個預測類的例項，每行代表一個實際類的例項
from sklearn.metrics import confusion_matrix
print confusion_matrix(classifier.predict(test), t_test)
"""
如果我們牢記所有正確的猜測都在表格的對角線上，那麼觀測表格的錯誤就很容易了，即
對角線以外的非零值
"""

# 可以展示分類器效能的完整報告的方法也是很好用的
from sklearn.metrics import classification_report
print classification_report(classifier.predict(test), t_test, target_names=["setosa", "versicolor", "virginica"])

#     ------ 聚類 ------
from sklearn.cluster import KMeans
kmeans = KMeans(n_clusters=3, init="random")
kmeans.fit(data)

c = kmeans.predict(data)

# 估計群集的結果，與使用完整性得分和同質性得分計算而得的標籤作比較
from sklearn.metrics import completeness_score, homogeneity_score
print completeness_score(t, c)
print homogeneity_score(t, c)
"""
當大部分資料點屬於一個給定的類並且屬於同一個群集，那麼完整性得分就趨向於 1
當所有群集都幾乎只包含某個單一類的資料點時同質性得分就趨向於 1
"""

# 把叢集視覺化並和帶有真實標籤的做視覺化比較
figure()
subplot(211)
plot(data[t == 1, 0], data[t == 1, 2], "bo")
plot(data[t == 2, 0], data[t == 2, 2], "ro")
plot(data[t == 3, 0], data[t == 3, 2], "go")
subplot(212)
plot(data[c == 1, 0], data[c == 1, 2], "bo", alpha=.7)
plot(data[c == 2, 0], data[c == 2, 2], "ro", alpha=.7)
plot(data[c == 0, 0], data[c == 0, 2], "go", alpha=.7)
# show()
close()

#     ------ 迴歸 ------
# 為了應用線性迴歸，我們需要建立一個由上所述的綜合資料集
from numpy.random import rand
x = rand(40, 1)
y = x*x*x + rand(40, 1) / 5

from sklearn.linear_model import LinearRegression
linreg = LinearRegression()
linreg.fit(x, y)

# 我們可以通過把擬合線和實際資料點畫在同一幅圖上來評估結果
from numpy import linspace, matrix
xx = linspace(0, 1, 40)
plot(x, y, "o", xx, linreg.predict(matrix(xx).T), "--r")
# show()
close()

# 我們還可以使用均方誤差來量化模型和原始資料的擬合度
from sklearn.metrics import mean_squared_error
print mean_squared_error(linreg.predict(x), y)
"""
該指標度量了預期的擬合線和真實資料之間的距離平方，當擬合線很完美時該值為 0
"""

#     ------ 相關 ------
"""
我們通過研究相關性來理解成對的變數之間是否相關，相關性的強弱。此類分析幫助我們精確定位被依賴的重要變數。
最好的相關方法是皮爾遜積矩相關係數，它是由兩個變數的協方差除以他們的標準差的乘機計算而來
"""

from numpy import corrcoef
corr = corrcoef(data.T)
print corr

from pylab import pcolor, colorbar, xticks, yticks, close
from numpy import arange
pcolor(corr)
colorbar()
xticks(arange(0.5, 4.5), ['sepal length',  'sepal width', 'petal length', 'petal width'],rotation=-20)
yticks(arange(0.5,4.5),['sepal length',  'sepal width', 'petal length', 'petal width'],rotation=-20)
# show()
close()

#    ------ 降維 ------
# 最著名的降維技術之一就是主成分分析
from sklearn.decomposition import PCA
pca = PCA(n_components=2)

pcad = pca.fit_transform(data)

plot(pcad[target == "setosa", 0], pcad[target == "setosa", 1], "bo")
plot(pcad[target == "versicolor", 0], pcad[target == "versicolor", 1], "ro")
plot(pcad[target == "virginica", 0], pcad[target == "virginica", 1], "go")
# show()
close()

# PCA 將空間資料方差最大化，我們可以通過方差比判斷 PCs 包含的資訊量
print pca.explained_variance_ratio_
# 輸出：[ 0.92461621  0.05301557]
# 現在我們知道第一個 PC 佔原始資料的 92% 的資訊量而第二個佔剩下的 5%，我們還可以輸出在轉化過程中
# 丟失的資訊量
print 1 - sum(pca.explained_variance_)
# 此時我們可以是應用逆變換還原原始資料
data_inv = pca.inverse_transform(pcad)
# 可以證明的是，由於資訊丟失逆變換不能給出準確的原始資料，我們可以估算逆變換的結果和原始資料的相似度
print abs(sum(sum(data - data_inv)))

# 通過改變主成分的數值來計算我們能夠覆蓋多少資訊量是很有趣的
for i in range(1, 5):
    pca = PCA(n_components=i)
    pca.fit(data)
    print sum(pca.explained_variance_ratio_) * 100, "%"

#     ------ 網路挖掘 ------
# 通常我們分析的資料是以網路結構儲存的，我們可以使用點和邊描述之間的關係
# 本章中我們將會介紹分析此類資料的基本步驟，稱為圖論，一個幫助我們創造、處理和研究網路的類庫
# 尤其我們將會介紹如何使用特定方法建立有意義的資料視覺化，以及如何建立一組關聯稠密的點
# 使用圖論可以讓我們很容易的匯入用於描述資料結構的最常用結構
import networkx as nx
G = nx.read_gml("lesmiserables.gml", relabel=True)      # networkx 必須要下載 1.9.1 版本才行
# 在上述程式碼我們匯入了《悲慘世界》同時出現的單片語成的網路，可以通過https://gephi.org/datasets/lesmiserables.gml.zip免費
# 下載，資料以GML格式儲存。我們還可以使用下面的命令匯入並可視化網路：
nx.draw(G, node_size=0, edge_color="b", alpha=.2, font_size=7)