包版本信息¶
In [3]:
!pip freeze | grep numpy
!pip freeze | grep pandas
!pip freeze | grep scikit-learn
数据预处理¶
In [57]:
import numpy as np
import pandas as pd
## 隐藏warning信息
import warnings
warnings.filterwarnings('ignore')
由于数据集以txt文件存储，使用Pandas的read_table()函数读入数据，分隔符为\t。
In [58]:
data = pd.read_table('seed.txt', sep='\t', names = ['x1', 'x2', 'x3', 'x4', 'x5', 'x6', 'x7','label'])
data.head()
Out[58]:
然后，我们将data['label']的取值减1, 标签取值范围为[0, 1, 2]。
In [59]:
data['label'] = data['label'] - 1
data.head()
Out[59]:
In [60]:
data.shape
Out[60]:
In [61]:
# 检查有无缺失值
np.sum(data.isnull())
Out[61]:
In [62]:
# 标签label的取值1，2，3
data['label'].value_counts()
Out[62]:
In [63]:
class_label = np.unique(data['label'])
class_label
Out[63]:
In [64]:
data.info()
In [65]:
data.describe()
Out[65]:
In [19]:
from sklearn.cross_validation import cross_val_score
from sklearn.model_selection import GridSearchCV
from sklearn.model_selection import train_test_split
from sklearn.linear_model.logistic import LogisticRegression
from sklearn.metrics import classification_report, accuracy_score, confusion_matrix
In [20]:
def evaluate(pred,test_y):
import seaborn as sns
import matplotlib.pyplot as plt
import sklearn.metrics as metrics
# 输出分类的准确率
print("Accuracy: %.4f" % (metrics.accuracy_score(test_y,pred)))
# 输出衡量分类效果的各项指标
print(classification_report(test_y, pred))
# 更直观的，我们通过seaborn画出混淆矩阵
%matplotlib inline
plt.figure(figsize=(6,4))
colorMetrics = metrics.confusion_matrix(test_y,pred)
# 坐标y代表test_y，即真实的类别，坐标x代表估计出的类别pred
sns.heatmap(colorMetrics, annot=True, fmt='d', xticklabels=[1,2,3],yticklabels=[1,2,3])
sns.plt.show()
In [23]:
## 数据准备
X = data.iloc[:, :-1]
y = data['label']
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.25)
In [24]:
lr = LogisticRegression()
param_dict = {'C': [0.001, 0.01, 0.1, 1, 10]}
gs = GridSearchCV(lr, param_grid = param_dict, scoring='accuracy')
gs.fit(X_train, y_train)
Out[24]:
In [25]:
pd.DataFrame(gs.cv_results_)
Out[25]:
In [26]:
print gs.best_params_
print gs.best_score_
In [27]:
lr_cls = gs.best_estimator_
score = cross_val_score(lr_cls, X_train, y_train, cv = 10, scoring = 'accuracy')
print np.mean(score)
In [28]:
evaluate(y_test, lr_cls.predict(X_test))
- 计算Expectation步
在这一步，我们计算每条数据记录$x_{n}$分属K个标签的概率$\gamma_{n}(z_{k}) = p(z_{k} = 1 \vert x_{n}, \theta)$。
In [33]:
'''expectation_step函数
输入参数：X: 数据集
pis: p(z = 1)的概率
mus: 高斯分布的均值
sigmas: 高斯分布的协方差
返回参数： $x_{n}$分属K个标签的概率矩阵
'''
def expectation_step(X, pis, mus, sigmas):
import numpy as np
from scipy.stats import multivariate_normal
gammas = np.array([value[0]*multivariate_normal(mean = value[1], cov = value[2]).pdf(X) for value in zip(pis, mus, sigmas)])
return gammas / gammas.sum(axis = 0)
- 计算Maximization步
在这一步，我们以完全数据的对数极大似然函数为目标，更新参数pis, mus, sigmas
In [34]:
'''maximization_step函数
输入参数：X: 数据集
gammas：$x_{n}$分属K个标签的概率矩阵
pis: p(z = 1)的概率
mus: 高斯分布的均值
sigmas: 高斯分布的协方差
'''
def maximization_step(X, gammas, pis, mus, sigmas):
# 更新pis
pis = gammas.sum(axis = 1) / gammas.shape[1]
# 更新mus
mus = np.multiply(1.0/gammas.sum(axis = 1).reshape(-1, 1), np.dot(gammas, X))
# 更新sigmas
for k in range(len(sigmas)):
centered_X = X - mus[k]
sigmas[k] = np.dot(np.multiply(gammas[k], centered_X.T), centered_X) / gammas[k].sum()
更新最大似然函数的下界lower_bound。
In [35]:
'''cal_lower_bound函数
输入参数：X: 数据集
gammas：$x_{n}$分属K个标签的概率矩阵
pis: p(z = 1)的概率
mus: 训练得到的高斯分布的均值
sigmas: 训练得到的高斯分布的协方差
返回参数：对数最大似然函数的下界
'''
def cal_lower_bound(X, gammas, pis, mus, sigmas):
from scipy.stats import multivariate_normal
lower_bound = 0
for gamma, pi, mu, sigma in zip(gammas, pis, mus, sigmas):
lower_bound += np.multiply(gamma, np.log(pi*multivariate_normal(mu, sigma).pdf(X)))
return lower_bound.sum()
接下来，我们实现train_GMM函数，计算各个高斯分布的统计特征。
In [36]:
'''train_GMM函数
输入参数：X: 数据集
k: 高斯分布的个数
max_iter: 迭代的最高次数
alpha: 最大似然函数下界收敛的阈值
返回参数：(mus, sigmas)以元组形式返回参数
mus: 高斯分布的均值 K times d
sigmas：高斯分布的协方差 d times d
'''
def train_GMM(X, K, max_iter, alpha):
indicator = 0
lower_bound = 0
# 初始化
N, d = X.shape
d = d - 1
# 初始化pi，一维数组，即 K times 1的向量
pis = np.random.random(K)
pis /= pis.sum()
# 初始化mu，二维数组，即 K times d的矩阵
#mus = np.random.random((K, d))
mus = np.array([X[X['label'] == num].mean()[:-1] for num in range(K)])
# 初始化sigmas 三维数组， 即长度为K，元素为 d times d 矩阵的数组
sigmas = np.array([np.eye(d)] * K)
X_hat = X.iloc[:,:-1].values
while indicator <= max_iter:
lower_bound_old = lower_bound
# E-step
gammas = expectation_step(X_hat, pis, mus, sigmas)
# M-step
maximization_step(X_hat, gammas, pis, mus, sigmas)
# Update 下界
lower_bound = cal_lower_bound(X_hat, gammas, pis, mus, sigmas)
# 检查终止条件
if abs(lower_bound - lower_bound_old) < alpha:
print 'reach stopping criterion!'
break
# 更新迭代次数
indicator += 1
return (mus, sigmas)
计算高斯分布的参数之后，我们需要判定样本来自于各个分布的概率大小。也就是说将样本代入分布的密度函数，计算其密度，取最大值所在的分布作为该样本的总体分布，并将样本归属于样本所在的类。
In [37]:
'''predict函数
输入参数： X：数据集
mus: 各个高斯分布的均值 K time d
sigmas: 各个高斯分布的协方差 K times 1 元素为d times d 矩阵
输出参数：数据集的各个数据记录所对应的标签(如0,1,2)
'''
def predict(X, mus, sigmas):
from scipy.stats import multivariate_normal
labels = np.zeros(X.shape[0])
for index, item in enumerate(X.iloc[:,:-1].values):
probs = [multivariate_normal(mu, sigma).pdf(item) for mu, sigma in zip(mus, sigmas)]
labels[index] = probs.index(max(probs))
return labels
算法验证¶
我们人为生成300个分别来自3个不同高斯分布的数据样本，组成验证数据集。
In [38]:
# 生成来自高斯混合分布的数据
# 数据集来自3个高斯分布，每个分布产生100个数据样本，数据集的维度为2
#np.random.seed(5)
mus = [[0, 4], [-2, 0], [4, 4]]
sigs = [[[3, 0], [0, 0.5]], [[1, 0], [0, 2]], [[1, 0], [0,1]]]
# 数据集X
multi_data = [np.random.multivariate_normal(mu, sig, 100) for mu, sig in zip(mus, sigs)]
data = np.vstack(multi_data)
## 数据集的类别标签为0，1，2
labels = [[value]*100 for value in range(3)]
## 把嵌套数组labels拉平
multi_labels = [y for x in labels for y in x]
# 数据集标签合并
final_data = []
for v1, v2 in zip(data, multi_labels):
list_v1 = v1.tolist()
list_v1.append(v2)
final_data.append(list_v1)
##np.random.shuffle(final_data)
把数据集转化为Pandas中的DataFrame格式。
In [39]:
final_data = pd.DataFrame(final_data, columns = ['x1', 'x2', 'label'])
final_data.head()
Out[39]:
测试数据集准备好之后，我们把final_data代入函数中进行计算。
In [40]:
mus, sigmas = train_GMM(final_data, 3, 1000, 1e-5)
In [41]:
# 根据训练集生成的均值
mus
Out[41]:
In [42]:
# 根据训练集生成的协方差
sigmas
Out[42]:
In [43]:
# 在训练集上预测的标签
label_pred = predict(final_data, mus, sigmas)
print label_pred
In [44]:
## 模型评价
evaluate(label_pred, final_data['label'].values)
使用sklearn的Gaussian函数¶
In [45]:
from sklearn.mixture import GaussianMixture
In [46]:
gmm_model = GaussianMixture(n_components = 3, covariance_type = 'diag', init_params='random', max_iter=200)
gmm_model.means_init = np.array([final_data[final_data['label'] == num].mean(axis = 0)[:-1] for num in range(3)])
gmm_model.fit(final_data.iloc[:,:-1])
y_pred = gmm_model.predict(final_data.iloc[:,:-1])
diff = y_pred - final_data['label'].values
evaluate(y_pred, final_data['label'].values)
接下来，我们查看训练得到的均值和协方差矩阵。
In [47]:
gmm_model.means_
Out[47]:
In [48]:
gmm_model.covariances_
Out[48]:
真实数据集¶
使用sklearn中的GaussianMixture模型，查看数据集是否适合使用高斯混合模型。
In [66]:
gmm_model = GaussianMixture(n_components = 3, covariance_type = 'diag', init_params='random', max_iter=200)
gmm_model.means_init = np.array([data[data['label'] == num].mean(axis = 0)[:-1] for num in range(3)])
gmm_model.fit(data.iloc[:,:-1])
y_pred = gmm_model.predict(data.iloc[:,:-1])
diff = y_pred - data['label'].values
evaluate(y_pred, data['label'].values)
准确率接近90%，数据集可以使用高斯混合模型来进行分类，也就是认为数据集的各个部分服从高斯分布。接下来，我们使用自定义的模型，去拟合数据集。
首先，使用数据集训练模型。我们将模型中高斯分布的个数设为3，迭代次数设为2000，终止条件的阈值设为1e-5。
In [71]:
mus, sigmas = train_GMM(data, 3, 2000, 1e-5)
train_GMM()函数返回各个高斯分布component的参数值。
In [74]:
## 均值
mus
Out[74]:
In [75]:
## 协方差矩阵
sigmas
Out[75]:
预测数据集的标签。
In [76]:
label_pred = predict(data, mus, sigmas)
label_pred
Out[76]:
In [73]:
evaluate(label_pred, data['label'].values)
