K-means聚类简介

前面介绍了不少监督学习的各类算法，那么这次我们就来认识和了解一下无监督学习。首先，K-means聚类算法就是一种无监督学习的算法。这就意味着他可以通过没有分类标记的数据学习出训练数据本身的分组信息。（其实他的本质就是自动的将属性/特征类似的数据归到统一类别）而其最后分成的组数则由参数K来决定。

算法基本步骤

1. 随机初始化K个簇心。用作标识各组数据。
2. 根据各个数据的特征/属性，将他们聚类到距离最接近的一个簇，并打上分类标识。
3. 根据加入的数据重新计算每个簇的质心（均值，Means）
4. 重复2，3步，直至最后所有数据被划分成k个簇。

从以后的步骤我们也可以充分地看出，再聚类完毕之前数据是没有分组和标识的，而是通过聚类来查找和分析被划分在一起的数据。

同时，每个簇的质心是定义组的特征的代表，我们往往可以通过一个簇的质心看出该分组的特征以及属性。从而判断出这个簇代表了一个怎么样的组。

应用领域

• 网站，App的用户群体聚类
• 信用卡诈骗监测（这个是利用了聚类算法的反原理，寻找离群值）
• 文章自动归类
• 图像归类
• 人群健康状态监测（同信用卡诈骗检测同原理）

此外，我们还可以通过监控数据点根据时间的变化而造成的组间切换现象找到状态变化的方法。

算法主要原理

（1）在数据划分方面，我们主要依赖的思路是：欧氏距离。即，每个簇的质心定义了一种聚类，然后每个点根据其到每个簇的质心的欧几里德距离的平方的大小，将其分配到最近的簇中。我们假设$c_i$是簇集合C中质心i，那么我们可以将其聚类的过程表达成下列形式：

$argmin_{c_i \in C} dist(c_i,x)^2$

（2）质心的更新，在这一步骤中，各个簇的质心将要被重新计算，通过加入新的数据再重新计算包括i新数据在内的的整个簇的均值以此作为新的质心。该过程可以表达为一下形式：

$c_i=\frac{1}{|S_i|}\sum_{x_i \in S_i }x_i$

在多次迭代之后，该算法可以保证收敛于一个结果，虽然有可能只是一个局部的最优解而已。为了使得最后结果更加稳定和可靠，我们也可以通过多次运行算法然后再在多个结果中取得一个均衡（为什么说多次计算可以得到更好的结果呢？ 因为每次随机初始化的质心都是不同的，这样使得每次的运行结果都不尽相同，为了得到一个更加稳定的结果和更好的鲁棒性，多次运行往往是一个简单但却有效的方法。）

有关K值的选择

在K-means聚类算法中，K值的选择往往是一个较为难以抉择的问题。最后的情况是，在你应用这个算法之前，已经知道数据应该分为多少个类别合适，或者将这种方法应用在一个半监督学习的场合，使得最后的分类结果有据可依。

那么，在一般我们不太清楚分类的情况之下，我们能做的就是通过为K值制定一个范围，然后在这个范围之上，我们尝试运行我们的分类器，然后根据分类效果来评判我们的k的取值。

另外一方面，聚类算法通常会不断降低各个样本点到质心之间的距离，随着k值的增加，“距离”更是一直下降。但是，当k值增大到一定程度时，虽然“距离”仍在下降，但下降的速率已经变得十分缓慢了。我们常常称这个点的k值为肘点。往往到达肘点的k值已经能够很好地将各个类别之间的差异信息给描述出来，所以我们可以在不清楚数据具体分类信息的情况下，选择肘点值来作为我们的k值。

Example

在这里我们通过Python的Sklearn库来看一下K-means的应用实例

Delivery Fleet Data 数据下载

数据读取

1
2
3
4
5
6
7
8
9

import pandas as pd
import matplotlib.pyplot as plt
df = pd.read_csv("data_1024.csv",sep='\t')
plt.scatter(df['Distance_Feature'],df['Speeding_Feature'])
plt.xlabel("Distance_Feature")
plt.ylabel("Speeding_Feature")
plt.title("Delivery Fleet Data")
plt.show()
df.head()

输出结果
Driver_ID Distance_Feature Speeding_Feature
0 3423311935 71.24 28.0
1 3423313212 52.53 25.0
2 3423313724 64.54 27.0
3 3423311373 55.69 22.0
4 3423310999 54.58 25.0

算法使用

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

import numpy as np
from sklearn.cluster import KMeans
X_train = df[['Distance_Feature','Speeding_Feature']]
kmeas = KMeans(n_clusters=2).fit(X_train)

color = [str(item/255.) for item in kmeas.labels_]
plt.scatter(df['Distance_Feature'],df['Speeding_Feature'],c=color)
plt.xlabel("Distance_Feature")
plt.ylabel("Speeding_Feature")
plt.title("Delivery Fleet Data")
plt.show()
kmeas = KMeans(n_clusters=4).fit(X_train)

color = [str(item/255.) for item in kmeas.labels_]
plt.scatter(df['Distance_Feature'],df['Speeding_Feature'],c=color)
plt.xlabel("Distance_Feature")
plt.ylabel("Speeding_Feature")
plt.title("Delivery Fleet Data")
plt.show()

输出结果

上面的结果展示了，K分别等于2和4时对数据的聚类的结果，由于这个算法使用还算简单，在这里我也不做赘述了。那么有关K-means的简介到这里也就结束了，谢谢大家的阅读。

机器学习之梯度下降

梯度下降(gradient descent)

是线性回归的一种(Linear Regression)，在Andrew Ng的machine learning前期课程中有着详细的讲解，首先给出一个关于课程中经典的例子:预测房屋的价格。

在上面的表格中，房间的个数以及房屋的面积都是输入变量而，价格就是我们要预测输出的值。其中，面积和房屋个数分别表示一个特征，我们在这将他们一起用X来表示。价格用Y来表示。同时表格的每一行表示成一个样本，具体表示为：$f: X_1=[24104,3] \to Y_1=400$。

那么我们具体的任务就可以概括成给定一个训练集合，学习出一个函数$f$，使得$f(x)$能够符合结果Y，或者说能够使得$f(x)$与Y的之间的误差最小。

我们可以用以下的式子来表达上述的关系：

$f(x) = \theta_1 x_1 + \theta_2 x_2 + \beta$

其中，$\theta 表示X映射成Y的权重，而x表示为各组特征。$
我们再分别用$x^{j},y^{j}来表示第J个样本$。这样我们就可以写出代价函数：

$J(\theta) = \frac {1}{2}\sum ^m_{i=0}(f_\theta(x^{i})-y^{(i)})^2$

然后我们要使得计算的结果无限接近真实值y的话，我们就要令得我们的代价函数（俗称的误差）最小化。要使得$J(\theta)$最小，即对$J(\theta)$进行求导操作，并使得其结果为0。

$\theta_j = \theta_j - \frac{\partial J(\theta)}{\partial \theta_j}$

当只有单个特征时，上式中的j表示系数（权重）的编号，右边的值赋给左边的变量后完成一次迭代过程。

而多个特征时，其迭代过程如下：

上式就是批梯度下降算法(batch gradient descent)，当上式收敛时则退出迭代，何为收敛，即前后两次迭代的值不再发生变化了。一般情况下，会设置一个具体的参数，当前后两次迭代差值小于该参数时候结束迭代。注意以下几点：

(1) a 即学习率（learning rate），决定的下降步伐，如果太小，则找到函数最小值的速度就很慢，如果太大，则可能会出现无法逼近最小值的现象；

(2) 初始点不同，获得的最小值也不同，因此梯度下降求得的只是局部最小值；

(3) 越接近最小值时，下降速度越慢；

(4) 计算批梯度下降算法时候，计算每一个θ值都需要遍历计算所有样本，当数据量的时候这是比较费时的计算。

批梯度下降算法的步骤可以归纳为以下几步：

(1)先确定向下一步的步伐大小，我们称为Learning rate ；

(2)任意给定一个初始值：θ向量，一般为0向量

(3)确定一个向下的方向，并向下走预先规定的步伐，并更新θ向量

(4)当下降的高度小于某个定义的值，则停止下降；

随机梯度下降算法

因为每次计算梯度都需要遍历所有的样本点。这是因为梯度是J(θ)的导数，而J(θ)是需要考虑所有样本的误差和 ，这个方法问题就是，扩展性问题，当样本点很大的时候，基本就没法算了。

所以接下来又提出了随机梯度下降算法(stochastic gradient descent )。随机梯度下降算法，每次迭代只是考虑让该样本点的J(θ)趋向最小，而不管其他的样本点，这样算法会很快，但是收敛的过程会比较曲折，整体效果上，大多数时候它只能接近局部最优解，而无法真正达到局部最优解。所以适合用于较大训练集的例子。

其具体流程用公式表达如下：

代码实现

批梯度下降算法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

#Training data set
#each element in x represents (x0,x1,x2)
x = [(1,0.,3) , (1,1.,3) ,(1,2.,3), (1,3.,2) , (1,4.,4)]
#y[i] is the output of y = theta0 * x[0] + theta1 * x[1] +theta2 * x[2]
y = [95.364,97.217205,75.195834,60.105519,49.342380]


epsilon = 0.000001
#learning rate
alpha = 0.001
diff = [0,0]
error1 = 0
error0 =0
m = len(x)

#init the parameters to zero
theta0 = 0
theta1 = 0
theta2 = 0
sum0 = 0
sum1 = 0
sum2 = 0
while True:

    #calculate the parameters
    for i in range(m):
        #begin batch gradient descent
        diff[0] = y[i]-( theta0 + theta1 * x[i][1] + theta2 * x[i][2] )
        sum0 = sum0 + alpha * diff[0]* x[i][0]
        sum1 = sum1 + alpha * diff[0]* x[i][1]
        sum2 = sum2 + alpha * diff[0]* x[i][2]
        #end  batch gradient descent
    theta0 = sum0;
    theta1 = sum1;
    theta2 = sum2;
    #calculate the cost function
    error1 = 0
    for lp in range(len(x)):
        error1 += ( y[i]-( theta0 + theta1 * x[i][1] + theta2 * x[i][2] ) )**2/2

    if abs(error1-error0) < epsilon:
        break
    else:
        error0 = error1

    print(' theta0 : %f, theta1 : %f, theta2 : %f, error1 : %f'%(theta0,theta1,theta2,error1))

print ('Done: theta0 : %f, theta1 : %f, theta2 : %f'%(theta0,theta1,theta2))

输出结果：
theta0 : 0.377225, theta1 : 0.625295, theta2 : 1.120912, error1 : 4405.869965
theta0 : 0.729497, theta1 : 1.193311, theta2 : 2.164098, error1 : 3094.652486
theta0 : 1.058680, theta1 : 1.708424, theta2 : 3.135362, error1 : 2089.261446
theta0 : 1.366497, theta1 : 2.174683, theta2 : 4.040070, error1 : 1335.974025
theta0 : 1.654541, theta1 : 2.595831, theta2 : 4.883186, error1 : 789.589683
theta0 : 1.924288, theta1 : 2.975327, theta2 : 5.669299, error1 : 412.137681
theta0 : 2.177098, theta1 : 3.316371, theta2 : 6.402654, error1 : 171.776368
theta0 : 2.414234, theta1 : 3.621921, theta2 : 7.087177, error1 : 41.856096
theta0 : 2.636861, theta1 : 3.894714, theta2 : 7.726499, error1 : 0.121739
theta0 : 2.846057, theta1 : 4.137277, theta2 : 8.323976, error1 : 28.034282
theta0 : 3.042819, theta1 : 4.351950, theta2 : 8.882714, error1 : 110.193963
……
theta0 : 98.101057, theta1 : -13.028753, theta2 : 1.133773, error1 : 3.473679
theta0 : 98.101058, theta1 : -13.028753, theta2 : 1.133773, error1 : 3.473678
theta0 : 98.101058, theta1 : -13.028753, theta2 : 1.133773, error1 : 3.473677
theta0 : 98.101059, theta1 : -13.028753, theta2 : 1.133773, error1 : 3.473676
theta0 : 98.101059, theta1 : -13.028753, theta2 : 1.133773, error1 : 3.473675
theta0 : 98.101060, theta1 : -13.028753, theta2 : 1.133772, error1 : 3.473674
theta0 : 98.101061, theta1 : -13.028753, theta2 : 1.133772, error1 : 3.473673
Done: theta0 : 98.101061, theta1 : -13.028753, theta2 : 1.133772

由上面的输出结果可以知道，梯度下降算法会在迭代一定次数后收敛于一个较少的损失值（即error）然而这并不是最优解而只是一个局部最小值（因为我们也可以从结果看到 theta0 : 2.636861, theta1 : 3.894714, theta2 : 7.726499, error1 : 0.121739，这项数据的最终error反而优于我们的最终解）。

随机梯度下降算法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45

#Training data set
#each element in x represents (x0,x1,x2)
x = [(1,0.,3) , (1,1.,3) ,(1,2.,3), (1,3.,2) , (1,4.,4)]
#y[i] is the output of y = theta0 * x[0] + theta1 * x[1] +theta2 * x[2]
y = [95.364,97.217205,75.195834,60.105519,49.342380]


epsilon = 0.0001
#learning rate
alpha = 0.01
diff = [0,0]
error1 = 0
error0 =0
m = len(x)


#init the parameters to zero
theta0 = 0
theta1 = 0
theta2 = 0

while True:

    #calculate the parameters
    for i in range(m):

        diff[0] = y[i]-( theta0 + theta1 * x[i][1] + theta2 * x[i][2] )

        theta0 = theta0 + alpha * diff[0]* x[i][0]
        theta1 = theta1 + alpha * diff[0]* x[i][1]
        theta2 = theta2 + alpha * diff[0]* x[i][2]

    #calculate the cost function
    error1 = 0
    for lp in range(len(x)):
        error1 += ( y[i]-( theta0 + theta1 * x[i][1] + theta2 * x[i][2] ) )**2/2

    if abs(error1-error0) < epsilon:
        break
    else:
        error0 = error1

    print (' theta0 : %f, theta1 : %f, theta2 : %f, error1 : %f'%(theta0,theta1,theta2,error1))

print ('Done: theta0 : %f, theta1 : %f, theta2 : %f'%(theta0,theta1,theta2))

输出结果：
theta0 : 2.782632, theta1 : 3.207850, theta2 : 7.998823, error1 : 7.508687
theta0 : 4.254302, theta1 : 3.809652, theta2 : 11.972218, error1 : 813.550287
theta0 : 5.154766, theta1 : 3.351648, theta2 : 14.188535, error1 : 1686.507256
theta0 : 5.800348, theta1 : 2.489862, theta2 : 15.617995, error1 : 2086.492788
theta0 : 6.326710, theta1 : 1.500854, theta2 : 16.676947, error1 : 2204.562407
theta0 : 6.792409, theta1 : 0.499552, theta2 : 17.545335, error1 : 2194.779569
theta0 : 7.223066, theta1 : -0.467855, theta2 : 18.302105, error1 : 2134.182794
theta0 : 7.630213, theta1 : -1.384304, theta2 : 18.982980, error1 : 2056.719790
……
theta0 : 97.986505, theta1 : -13.221170, theta2 : 1.257223, error1 : 1.553680
theta0 : 97.986620, theta1 : -13.221169, theta2 : 1.257186, error1 : 1.553579
theta0 : 97.986735, theta1 : -13.221167, theta2 : 1.257150, error1 : 1.553479
theta0 : 97.986849, theta1 : -13.221166, theta2 : 1.257113, error1 : 1.553379
theta0 : 97.986963, theta1 : -13.221165, theta2 : 1.257077, error1 : 1.553278
Done: theta0 : 97.987078, theta1 : -13.221163, theta2 : 1.257041

通过上述批梯度下降和随机梯度下降算法代码的对比，不难发现两者的区别：
随机梯度下降算法在迭代的时候，每迭代一个新的样本，就会更新一次所有的theta参数。
因此当样本数量很大时候，批梯度得做完所有样本的计算才能更新一次theta，从而花费的时间远大于随机梯度下降。但是随机梯度下降过早的结束了迭代，使得它获取的值只是接近局部最优解，而并非像批梯度下降算法那样就是局部最优解。

Python 实现推荐系统

引言

最早的推荐系统应该是亚马逊为了提升长尾货物的用户抵达率而发明的。已经有数据证明，长尾商品的销售额以及利润总和与热门商品是基本持平的。亚马逊网站上在线销售的商品何止百万，但首页能够展示的商品数量又极其有限，给用户推荐他们可能喜欢的商品就成了一件非常重要的事情。当然，商品搜索也是一块大蛋糕，亚马逊的商品搜索早已经开始侵蚀谷歌的核心业务了。

从这些例子之中，我们可以看到我们能够使用许多不同的方式来搜集兴趣偏好。有时候，这些数据可能来自人们购买的商品，以及这些商品关联的评价信息。我们可以利用一组算法从中挖掘，建立几个有意思的推荐系统。

推荐系统的简介

两种最普遍的推荐系统的类型是基于内容的推荐系统和协同过滤推荐系统（CF）。协同过滤基于用户对产品的态度产生推荐，基于内容的推荐系统基于物品属性的相似性进行推荐。CF可以分为基于内存的协同过滤和基于模型的协同过滤。

基于内容推荐

基于内容的推荐（Content-based Recommendation），它是建立在项目的内容信息上作出推荐的，而不需要依据用户对项目的评价意见，更多地需要用机器学习的方法从关于内容的特征描述的事例中得到用户的兴趣资料。在基于内容的推荐系统中，项目或对象是通过相关的特征的属性来定义，系统基于用户评价对象的特征，学习用户的兴趣，考察用户资料与待预测项目的相匹配程度。用户的资料模型取决于所用学习方法，常用的有决策树、神经网络和基于向量的表示方法等。

基于内容推荐方法的优点是：

• 不需要其它用户的数据，没有冷开始问题和稀疏问题。
• 能为具有特殊兴趣爱好的用户进行推荐。
• 能推荐新的或不是很流行的项目，没有新项目问题。
• 通过列出推荐项目的内容特征，可以解释为什么推荐那些项目。
• 已有比较好的技术，如关于分类学习方面的技术已相当成熟。

缺点是:

要求内容能容易抽取成有意义的特征，要求特征内容有良好的结构性，并且用户的口味必须能够用内容特征形式来表达，不能显式地得到其它用户的判断情况。

协同过滤推荐

协同过滤推荐（Collaborative Filtering Recommendation），是推荐系统中应用最早和最为成功的技术之一。它一般采用最近邻技术，利用用户的历史喜好信息计算用户之间的距离，然后利用目标用户的最近邻居用户对商品评价的加权评价值来预测目标用户对特定商品的喜好程度，系统从而根据这一喜好程度来对目标用户进行推荐。协同过滤最大优点是对推荐对象没有特殊的要求，能处理非结构化的复杂对象，如音乐、电影。

协同过滤是基于这样的假设：为一用户找到他真正感兴趣的内容的好方法是首先找到与此用户有相似兴趣的其他用户，然后将他们感兴趣的内容推荐给此用 户。其基本思想非常易于理解，在日常生活中，我们往往会利用好朋友的推荐来进行一些选择。协同过滤正是把这一思想运用到电子商务推荐系统中来，基于其他用 户对某一内容的评价来向目标用户进行推荐。

基于协同过滤的推荐系统可以说是从用户的角度来进行相应推荐的，而且是自动的即用户获得的推荐是系统从购买模式或浏览行为等隐式获得的，不需要用户努力地找到适合自己兴趣的推荐信息，如填写一些调查表格等。

和基于内容的过滤方法相比，

协同过滤具有如下的优点：

• 能够过滤难以进行机器自动内容分析的信息，如艺术品，音乐等。
• 共享其他人的经验，避免了内容分析的不完全和不精确，并且能够基于一些复杂的，难以表述的概念（如信息质量、个人品味）进行过滤。
• 有推荐新信息的能力。可以发现内容上完全不相似的信息，用户对推荐信息的内容事先是预料不到的。这也是协同过滤和基于内容的过滤一个较大的差别，基于内容的过滤推荐很多都是用户本来就熟悉的内容，而协同过滤可以发现用户潜在的但自己尚未发现的兴趣偏好。
• 能够有效的使用其他相似用户的反馈信息，较少用户的反馈量，加快个性化学习的速度。

缺点是：

最典型的问题有，
稀疏问题（Sparsity）
可扩展问题（Scalability）
冷启动问题

实例应用

我们下面通过使用MovieLens数据集（这是一个有关电影评分的经典数据集， 其中Movielens-100k和movielens-1M有用户对电影的打分，电影的title、genre、IMDB链接、用户的gender、age、occupation、zip code。movielens-10M中还有用户对电影使用的tag信息。）

数据集的下载地址

数据读取

1
2
3
4
5
6
7
8
9

import numpy as np
import pandas as pd
header = ['user_id','item_id','rating','timestamp']
df = pd.read_csv('/home/ef/Desktop/ml-100k/u.data',
                sep='\t',names=header)
print(df.head())
n_users = df.user_id.unique().shape[0]
n_items = df.item_id.unique().shape[0]
print('Number of users = ' + str(n_users) + ' | Number of movies = ' + str(n_items))

输出结果：
user_id item_id rating timestamp
0 196 242 3 881250949
1 186 302 3 891717742
2 22 377 1 878887116
3 244 51 2 880606923
4 166 346 1 886397596
Number of users = 943 | Number of movies = 1682

数据划分

使用scikit-learn库来划分数据集

1
2

from sklearn import cross_validation as cv
X_train,x_test = cv.train_test_split(df,test_size=0.2)

基于内存的协同过滤

基于内存的协同过滤方法可以分为两个部分：

• 用户－产品协同过滤，即：基于用户的协同过滤推荐
• 产品－产品协同过滤，即：基于物品的协同过滤推荐

基于用户的协同过滤推荐，将选取一个特定的用户，基于打分的相似性发现类似于该用户的用户，并推荐那些相似用户喜欢的产品。基于物品的协同过滤推荐会选取一个产品，发现喜欢该产品的用户，并找到这些相似用户还喜欢的其它产品。

基于用户的协同过滤推荐：“喜欢这东西的人也喜欢……”
基于物品的协同过滤推荐：“像你一样的人也喜欢（买个该商品的还买了）……”

在这两种情况下，从整个数据集构建一个用户产品矩阵。

用户产品矩阵的例子：

计算相似性，并创建一个相似性矩阵。
通常用于推荐系统中的距离矩阵是余弦相似性，其中，打分被看成n维空间中的向量，而相似性是基于这些向量之间的角度进行计算的。用户a和b的余弦相似性可以用下面的公式进行计算：

$Similar_u^{cos}(U_a,U_b)=\frac{U_a*U_b}{||U_a|| * ||U_b||} = \frac {\sum x_{a,m}x_{b,m}}{\sqrt {\sum x^2_{a,m}\sum x^2_{b,m} }}$

同时，物品a和b的相似读也可以写成以上形式。

创建用户产品矩阵，针对测试数据和训练数据，创建两个矩阵：

1
2
3
4
5
6
7
8
9
10
11
12
13
14

from sklearn import cross_validation as cv
X_train,x_test = cv.train_test_split(df,test_size=0.2)

train_data_matrix = np.zeros((n_users,n_items))
for line in X_train.itertuples():
    train_data_matrix[line[1]-1, line[2]-1] = line[3]
test_data_matrix = np.zeros((n_users, n_items))
for line in x_test.itertuples():
    test_data_matrix[line[1]-1, line[2]-1] = line[3]

## 通过余弦相似度计算用户和物品的相似度
from sklearn.metrics.pairwise import pairwise_distances
user_similarity = pairwise_distances(train_data_matrix, metric = "cosine")
item_similarity = pairwise_distances(train_data_matrix.T, metric = "cosine")

评估

这里采用均方根误差(RMSE)来度量预测评分的准确性,可以使用sklearn的mean_square_error(MSE)函数，其中RMSE仅仅是MSE的平方根。

1
2
3
4
5
6
7
8
9
10
11

from sklearn.metrics import mean_squared_error
from math import sqrt
def rmse(prediction, ground_truth):
    prediction = prediction[ground_truth.nonzero()].flatten()
    ground_truth = ground_truth[ground_truth.nonzero()].flatten()
    return sqrt(mean_squared_error(prediction, ground_truth))

print ('User based CF RMSE: ' + str(rmse(user_prediction, test_data_matrix)))
print ('Item based CF RMSe: ' + str(rmse(item_prediction, test_data_matrix)))
print ('User based CF RMSE: ' + str(rmse(user_prediction, test_data_matrix)))
print ('Item based CF RMSe: ' + str(rmse(item_prediction, test_data_matrix)))

输出结果：
User based CF RMSE: 3.087357412872858
Item based CF RMSe: 3.437038163412728

但从结果来看，算法的推荐效果还算理想。

基于模型的协同过滤

基于模型的协同过滤是基于矩阵分解(MF)的，矩阵分解广泛应用于推荐系统中，它比基于内存的CF有更好的扩展性和稀疏性。MF的目标是从已知的评分中学习用户的潜在喜好和产品的潜在属性，随后通过用户和产品的潜在特征的点积来预测未知的评分。
这是一个典型的机器学习的问题，可以将已有的用户喜好信息作为训练样本，训练出一个预测用户喜好的模型，这样以后用户在进入系统，可以基于此模型计算推荐。这种方法的问题在于如何将用户实时或者近期的喜好信息反馈给训练好的模型，从而提高推荐的准确度。

SVD

SVD即：奇异值分解（Singular value decomposition）奇异值分解是线性代数中一种重要的矩阵分解，在信号处理、统计学等领域有重要应用。奇异值分解在某些方面与对称矩阵或Hermite矩阵基于特征向量的对角化类似。然而这两种矩阵分解尽管有其相关性，但还是有明显的不同。对称阵特征向量分解的基础是谱分析，而奇异值分解则是谱分析理论在任意矩阵上的推广。

在矩阵M的奇异值分解中

$M=UΣV^T$

• U的列(columns)组成一套对M的正交”输入”或”分析”的基向量。这些向量是MM*的特征向量。
• V的列(columns)组成一套对M的正交”输出”的基向量。这些向量是M*M的特征向量。
• Σ对角线上的元素是奇异值，可视为是在输入与输出间进行的标量的”膨胀控制”。这些是M*M及MM*的奇异值，并与U和V的列向量相对应。

那么矩阵M就可以分解成U，Σ，V。
U矩阵表示对应于隐藏特性空间中的用户的特性矩阵，而V矩阵表示对应于隐藏特性空间中的产品的特性矩阵。

现在，我们可以通过U，Σ和$V^T$的点积进行预测了。

1
2
3
4
5
6
7
8
9
10

# 计算矩阵的稀疏度
sparsity = round(1.0 - len(df) / float(n_users*n_items),3)
print('The sparsity level of MovieLen100K is ' + str(sparsity * 100) + '%')

import scipy.sparse as sp
from scipy.sparse.linalg import svds
u, s, vt = svds(train_data_matrix, k = 20)
s_diag_matrix = np.diag(s)
x_pred = np.dot(np.dot(u,s_diag_matrix),vt)
print('User-based CF MSE: ' + str(rmse(x_pred, test_data_matrix)))

输出结果：
The sparsity level of MovieLen100K is 93.7%
User-based CF MSE: 2.63901831155016

最后，我们可以发现，采用SVD分解矩阵的基于模型的协同过滤方法在推荐的表现上更胜一筹。

关联分析

相关概念

作为经典的机器学习算法来说，关联分析并不算复杂。简单来说，从大规模数据集中寻找物品间的隐含关系被称作关联分析(association analysis)或者关联规则学习(association rule learning)。而其中关联分析的关键就是在众多的物品组成的项集(itemset)中找到那些经常一起出现的物品集合，也就是我们所谓的频繁项集。由于寻找物品的不同组合是一项非常耗时的任务，所需的计算代价也非常的高，这样显然用蛮力的搜索方法是不可取的。为了解决这样的问题，我们需要用到Apriori算法来解决。

有关关联分析的一些应用

1. 购物篮分析，这是关联分析最经典的一项应用。我们耳熟能详的尿布与啤酒就是来源于此（虽然这个故事的真实性有待考据，但这并不影响它成为一个营销界的小故事）。
2. 公众热点发现
3. 新闻及流行趋势挖掘
4. 搜索引擎推荐
5. 发现毒蘑菇的相似特征（该例子取自《机器学习实战》一书）

Apriori算法

优点：易编码实现
缺点：在大数据集上计算效果较慢

既然关联分析是一种在大规模数据集中寻找有趣的关系的任务。这些关系可以有两种形式：一种就是我们上面所说的频繁项集(frequent itemsets);另外一种则是关联规则(association rules)。关联规则，其暗示了两种物品之间可能存在很强的关系。

那么我们该如何定义有趣的关系呢？以及，频繁项集中的频繁有怎么划定呢？虽然，许多概念都可以说明这些问题，但是最重要及最通用的莫过于其支持度和可信度。

支持度(support)：一个项集的支持度被定义为数据集包含该项集的记录比例。

举个例子:

如上图所示，因为5项记录中有4条包含了{豆奶}，所以其支持度为4/5;依此类推，{豆奶，尿布}的支持读为3/5。

可信度或置信度(confidence)：这个概念是针对一条{尿布} -> {葡萄酒}这样的关联规则来定义的。

$置信度({尿布} -> {葡萄酒} )= 支持度({尿布,葡萄酒}) / 支持度({尿布})$

支持度和可信度是用来量化关联分析是否成功的方法。

Apriori原理

如果计算一个集合中的频繁项集的支持度，首先需要遍历全部可能的项集，比如针对一个包含了4个产品的集合，那么购买该集合产品的可能集合数目为2^4-1=15，而针对包含N项的集合，就需要遍历2^N-1。显然，这样的计算量很大。为了降低所需计算时间，就需要我们所谓的Apriori原理。

其基本思路：如果某个项集是频繁的，那么它的所有子集也是频繁的。该定理的逆反定理为：如果某一个项集是非频繁的，那么它的所有超集（包含该集合的集合）也是非频繁的。Apriori原理的出现，可以在得知某些项集是非频繁之后，不需要计算该集合的超集，有效地避免项集数目的指数增长，从而在合理时间内计算出频繁项集。

算法基本流程

当集合中项的个数大于0时：
创建一个长度为k的候选项集列表
检查数据以确认每个项集都是频繁的
保留频繁项集，并构建k+1项组成的候选集列表

代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61

def load_dataset():
    "Load the sample dataset."
    return [[1, 3, 4], [2, 3, 5], [1, 2, 3, 5], [2, 5]]
def createC1(dataset):
    "Create a list of candidate item sets of size one."
    c1 = []
    for transaction in dataset:
        for item in transaction:
            if not [item] in c1:
                c1.append([item])
    c1.sort()
    #frozenset because it will be a ket of a dictionary.
    return list(map(frozenset, c1))
def scanD(dataset, candidates, min_support):
    "Returns all candidates that meets a minimum support level"
    sscnt = {}
    for tid in dataset:
        #print(tid)
        for can in candidates:
            if can.issubset(tid):
                sscnt.setdefault(can, 0)
                sscnt[can] += 1

    num_items = float(len(dataset))
    retlist = []
    support_data = {}
    for key in sscnt:
        support = sscnt[key] / num_items
        if support >= min_support:
            retlist.insert(0, key)
        support_data[key] = support
    return retlist, support_data
def aprioriGen(freq_sets, k):
    "Generate the joint transactions from candidate sets"
    retList = []
    lenLk = len(freq_sets)
    for i in range(lenLk):
        for j in range(i + 1, lenLk):
            L1 = list(freq_sets[i])[:k - 2]
            L2 = list(freq_sets[j])[:k - 2]
            L1.sort()
            L2.sort()
            if L1 == L2:
                retList.append(freq_sets[i] | freq_sets[j])
    return retList


def apriori(dataset, minsupport=0.5):
    "Generate a list of candidate item sets"
    C1 = createC1(dataset)
    D = list(map(set, dataset))
    L1, support_data = scanD(D, C1, minsupport)
    L = [L1]
    k = 2
    while (len(L[k - 2]) > 0):
        Ck = aprioriGen(L[k - 2], k)
        Lk, supK = scanD(D, Ck, minsupport)
        support_data.update(supK)
        L.append(Lk)
        k += 1
    return L, support_data

挖掘关联规则

上面，我们实现了用apriori算法找出了频繁项集。现在我们来尝试找出关联规则。那么什么是关联规则呢？我们可以举这样一个例子：假设我们有一个频繁项集是：{豆奶，莴苣}。那么则意味着一个人如果他已经购买了豆奶，那么他也有很大的几率会购买莴苣;但反过来却不一定成立：则如果一个人购买了莴苣，但还是不能说明他有很大几率会接着购买豆奶。

既然，在上面的例子中我们已经用支持度作为一种参考单位去量化一个集合的频繁关系，在这里我们同理也可以通过另一个概念来度量关联规则。没错，它就是我们已经在上面提过的可信度。根据可信度的计算法则和定义可知：

$置信度({尿布} -> {葡萄酒} )= 支持度({尿布,葡萄酒}) / 支持度({尿布})$

这也是我们上面已经提及过的问题。加上，在上面的例子之中我们已经得到了所有的频繁项集的支持度，那么为了求出各关联规则的可信度也就是剩下一个除法的问题而已。

如何从频繁项集中生成关联规则

一般对于一个频繁项集而言，如果我们直接生成关联规则的话，往往会有很多条。例如：我们对一个频繁项集{0，1，2，3}直接生成关联规则的话结果就如下图所示一样：

从上图我们就可以看出仅仅含4项的频繁项集就会生成出如此多的关联规则，那么如果频繁项集中含有10，20，甚至100个元素那么生成的候选关联规则的数据是相当庞大。

在此，我们可以参考apriori算法中挖掘频繁项集的思路。在发现频繁项集时我们通过最低支持度的方法来忽略掉一些项集支持度的计算;同理，在这里我们将这个思想应用在关联规则的生成之上，即：如果一条规则的不满足最低可信度的要求，那么其子集自然也不能满足最低可信度的要求。对于这部分规则的可信度计算我们可以直接忽略。如上图中的阴影部分。我们假设规则{0，1，2} -> 3不能满足最低可信度，那么，他的子集{1，2}->{0，3}，{0，2}->{1，3}等将不能满足最低可信度。

代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38

def generateRules(L, support_data, min_confidence=0.7):
    """Create the association rules
    L: list of frequent item sets
    support_data: support data for those itemsets
    min_confidence: minimum confidence threshold
    """
    rules = []
    for i in range(1, len(L)):
        for freqSet in L[i]:
            H1 = [frozenset([item]) for item in freqSet]
            print("freqSet", freqSet, 'H1', H1)
            if (i > 1):
                rules_from_conseq(freqSet, H1, support_data, rules, min_confidence)
            else:
                calc_confidence(freqSet, H1, support_data, rules, min_confidence)
    return rules


def calc_confidence(freqSet, H, support_data, rules, min_confidence=0.7):
    "Evaluate the rule generated"
    pruned_H = []
    for conseq in H:
        conf = support_data[freqSet] / support_data[freqSet - conseq]
        if conf >= min_confidence:
            print(freqSet - conseq, '--->', conseq, 'conf:', conf)
            rules.append((freqSet - conseq, conseq, conf))
            pruned_H.append(conseq)
    return pruned_H


def rules_from_conseq(freqSet, H, support_data, rules, min_confidence=0.7):
    "Generate a set of candidate rules"
    m = len(H[0])
    if (len(freqSet) > (m + 1)):
        Hmp1 = aprioriGen(H, m + 1)
        Hmp1 = calc_confidence(freqSet, Hmp1,  support_data, rules, min_confidence)
        if len(Hmp1) > 1:
            rules_from_conseq(freqSet, Hmp1, support_data, rules, min_confidence)

运行测试

1
2
3

dataset = load_dataset()
L,support_data = apriori(dataset)
generateRules(L,support_data)

输出结果
freqSet frozenset({3, 5}) H1 [frozenset({3}), frozenset({5})]
freqSet frozenset({1, 3}) H1 [frozenset({1}), frozenset({3})]
frozenset({1}) —> frozenset({3}) conf: 1.0
freqSet frozenset({2, 5}) H1 [frozenset({2}), frozenset({5})]
frozenset({5}) —> frozenset({2}) conf: 1.0
frozenset({2}) —> frozenset({5}) conf: 1.0
freqSet frozenset({2, 3}) H1 [frozenset({2}), frozenset({3})]
freqSet frozenset({2, 3, 5}) H1 [frozenset({2}), frozenset({3}), frozenset({5})]
Out[6]:
[(frozenset({1}), frozenset({3}), 1.0),
(frozenset({5}), frozenset({2}), 1.0),
(frozenset({2}), frozenset({5}), 1.0)]

从输出结果可以看出，其中{1} —> {3}，{5} —> {2}，{2} —> {5}这三条关联规则的可信度居然达到了1.0，而且{2}，{5}还存在一个相互关联的关系。除此之外，我们还可以发现：虽然{1} —> {3}这条关联的可信度达到了1.0，但是{3} —> {1}却不足0.7（因为我们的最低可信度标准为0.7，没有输出显示的规则，则意味着它们的可信度低于0.7)。为了可以看见输出更多的规则，我们可以降低下最低可信度的标准，再来执行一遍代码。

1
2
3

dataset = load_dataset()
L,support_data = apriori(dataset)
generateRules(L,support_data,0.5)

输出结果
freqSet frozenset({3, 5}) H1 [frozenset({3}), frozenset({5})]
frozenset({5}) —> frozenset({3}) conf: 0.6666666666666666
frozenset({3}) —> frozenset({5}) conf: 0.6666666666666666
freqSet frozenset({1, 3}) H1 [frozenset({1}), frozenset({3})]
frozenset({3}) —> frozenset({1}) conf: 0.6666666666666666
frozenset({1}) —> frozenset({3}) conf: 1.0
freqSet frozenset({2, 5}) H1 [frozenset({2}), frozenset({5})]
frozenset({5}) —> frozenset({2}) conf: 1.0
frozenset({2}) —> frozenset({5}) conf: 1.0
freqSet frozenset({2, 3}) H1 [frozenset({2}), frozenset({3})]
frozenset({3}) —> frozenset({2}) conf: 0.6666666666666666
frozenset({2}) —> frozenset({3}) conf: 0.6666666666666666
freqSet frozenset({2, 3, 5}) H1 [frozenset({2}), frozenset({3}), frozenset({5})]
frozenset({5}) —> frozenset({2, 3}) conf: 0.6666666666666666
frozenset({3}) —> frozenset({2, 5}) conf: 0.6666666666666666
frozenset({2}) —> frozenset({3, 5}) conf: 0.6666666666666666
Out[7]:
[(frozenset({5}), frozenset({3}), 0.6666666666666666),
(frozenset({3}), frozenset({5}), 0.6666666666666666),
(frozenset({3}), frozenset({1}), 0.6666666666666666),
(frozenset({1}), frozenset({3}), 1.0),
(frozenset({5}), frozenset({2}), 1.0),
(frozenset({2}), frozenset({5}), 1.0),
(frozenset({3}), frozenset({2}), 0.6666666666666666),
(frozenset({2}), frozenset({3}), 0.6666666666666666),
(frozenset({5}), frozenset({2, 3}), 0.6666666666666666),
(frozenset({3}), frozenset({2, 5}), 0.6666666666666666),
(frozenset({2}), frozenset({3, 5}), 0.6666666666666666)]

这次我们就可以在输出中看见{3} —> {1}了，它的可信度为0.6666666666666666。

到此，有关关联分析的apriori算法就告一段落了。后续还会有有关关联分析FP-growth算法的文章讲解。