月度归档：2018年08月

16.1 目的

回归分析（regression analysis)是确定两种或两种以上变量间相互依赖的定量关系的一种统计分析方法。运用十分广泛，回归分析按照涉及的变量的多少，分为一元回归和多元回归分析；按照因变量的多少，可分为简单回归分析和多重回归分析；按照自变量和因变量之间的关系类型，可分为线性回归分析和非线性回归分析。如果在回归分析中，只包括一个自变量和一个因变量，且二者的关系可用一条直线近似表示，这种回归分析称为一元线性回归分析。如果回归分析中包括两个或两个以上的自变量，且自变量之间存在线性相关，则称为多重线性回归分析.

16.2 模型

16.3 问题

假设我们有一组某企业广告费与销售额之间的对应关系的实际数据：对这些数据可视化结果如下：

可视化代码如下：

import matplotlib.pyplot as plt
import matplotlib.font_manager as fm
myfont = fm.FontProperties(fname='/home/hadoop/anaconda3/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/simhei.ttf')


x=[4,8,9,8,7,12,6,10,6,9]
y=[9,20,22,15,17,23,18,25,10,20]
#y1=[10,18,20,18,16,26,14,22,14,20]  #y=2x+2

plt.figure(figsize=(6,3))
plt.scatter(x,y)
#plt.plot(x,y1,color='blue',linestyle='dashed')
plt.xlabel("广告费(万)",fontproperties=myfont,size=12)
plt.ylabel("销售额（万）",fontproperties=myfont,size=12)
plt.show()

import matplotlib.pyplot as plt

import matplotlib.font_manager as fm

myfont = fm.FontProperties(fname='/home/hadoop/anaconda3/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/simhei.ttf')

x=[4,8,9,8,7,12,6,10,6,9]

y=[9,20,22,15,17,23,18,25,10,20]

#y1=[10,18,20,18,16,26,14,22,14,20] #y=2x+2

plt.figure(figsize=(6,3))

plt.scatter(x,y)

#plt.plot(x,y1,color='blue',linestyle='dashed')

plt.xlabel("广告费(万)",fontproperties=myfont,size=12)

plt.ylabel("销售额（万）",fontproperties=myfont,size=12)

plt.show()

如果该企业先增加广告费，如20万，那么它可能带来的销售额大概是多少？
解决这个问题的关键就是如何根据已有数据的总体趋势，拟合出一条直线，如图1.1中虚线所示，那么新给出的广告费就可预测其对应的销售额？如何拟合这条直线呢？

16.4 简单入手

假设我们要拟合的这条直线为一元一次直线，表达式为：
y=ax+b （1.1）
要拟合这条直线或确定这条直线，只要求出a和b即可。那如何得到a和b呢？
这条直线要满足什么条件才是最好的或最能体现这些样本的趋势？
如果能使预测值与实际值的差距最小，应该是一个不错的直线。
由此，想到最小二乘法，利用最小二乘法作为损失函数，然后通过使损失函数最小化来求参数a和b。

最损失函数的最小值，而且损失函数为凸函数，故可利用梯度为0来求出参数a和b
具体过程如下：

如此，a，b确定后，自然直线y=ax+b也就确定了，这样便可以根据新的x值，预测其y值了。

16.5、用Python解方程求出a和b

#根据y=ax+b随意写几个数
X=[4,8,9,8,7,12,6,10,6,9]
Y=[9,20,22,15,17,23,18,25,10,20]
#初始化以下序列
Xsum=0.0
X2sum=0.0
Ysum=0.0
XY=0.0
#得到序列元素个数或长度
n=len(X)
for i in range(n):
    Xsum+=X[i]
    Ysum+=Y[i]
    XY+=X[i]*Y[i]
    X2sum+=X[i]**2
a=(n*XY-Xsum*Ysum)/( n*X2sum-Xsum**2)
b=(Ysum-a*Xsum)/n
print('所求直线为: y=%f*x+%f' % (a,b) )

#根据y=ax+b随意写几个数

X=[4,8,9,8,7,12,6,10,6,9]

Y=[9,20,22,15,17,23,18,25,10,20]

#初始化以下序列

Xsum=0.0

X2sum=0.0

Ysum=0.0

XY=0.0

#得到序列元素个数或长度

n=len(X)

for i in range(n):

Xsum+=X[i]

Ysum+=Y[i]

XY+=X[i]*Y[i]

X2sum+=X[i]**2

a=(n*XY-Xsum*Ysum)/( n*X2sum-Xsum**2)

b=(Ysum-a*Xsum)/n

print('所求直线为: y=%f*x+%f' % (a,b) )

打印结果为：
所求直线为: y=1.980810*x+2.251599

16.6 用迭代方式求参数

以下利用迭代的方法求出参数a和b
直接通过解方程来求参数a和b，如果参数比较多，样本数也很多的情况下，计算量非常大，而且也不现实。因此，我们需要另辟蹊径！
还是以求函数y=x^2的最小值为例，我们可以通过其导数为0来求来求的是y最小的x值；我们也可以通过迭代的方式来求，先从某点开始，如x_0开始，然后沿梯度的方向，不断靠近最小值点，如下图所示：

为啥每次修改的值，都能往函数最小值那个方向前进呢？这里的奥秘在于，我们每次都是向函数的梯度的相反方向来修改。

什么是梯度呢？翻开大学高数课的课本，我们会发现梯度是一个向量，它指向函数值上升最快的方向。显然，梯度的反方向当然就是函数值下降最快的方向了。我们每次沿着梯度相反方向去修改的值，当然就能走到函数的最小值附近。

之所以是最小值附近而不是最小值那个点，是因为我们每次移动的步长不会那么恰到好处，有可能最后一次迭代走远了越过了最小值那个点。步长的选择是门手艺，如果选择小了，那么就会迭代很多轮才能走到最小值附近；如果选择大了，那可能就会越过最小值很远，收敛不到一个好的点上。

按照上面的讨论，我们就可以写出梯度下降算法的公式

X=[4,8,9,8,7,12,6,10,6,9]
Y=[9,20,22,15,17,23,18,25,10,20]

#初始化参数a、b
b=0.1
a=0.2
#alpha是学习率，代表的是迭代的步长
alpha=0.001


#定义线性模型函数
def h(i):
    return b+a*X[i]

def diff(i):
    return h(i)-Y[i]

#times 表示迭代1000次，那么基本上可以下降到梯度的最小值了
#每次迭代做的事情，就是式（1.13）、（1.14）两个公式，不停的去修改参数a，b
for times in range(1000):
    sum1=0
    sum2=0
    for i in range(m):
        sum1=sum1+diff(i)
        sum2=sum2+diff(i)*X[i]
        b=b-alpha*sum1
        a=a-alpha*sum2


print ("参数b的值:%.2f"%b)
print ("参数a的值:%.2f"%a)

X=[4,8,9,8,7,12,6,10,6,9]

Y=[9,20,22,15,17,23,18,25,10,20]

#初始化参数a、b

b=0.1

a=0.2

#alpha是学习率，代表的是迭代的步长

alpha=0.001

#定义线性模型函数

def h(i):

return b+a*X[i]

def diff(i):

return h(i)-Y[i]

#times 表示迭代1000次，那么基本上可以下降到梯度的最小值了

#每次迭代做的事情，就是式（1.13）、（1.14）两个公式，不停的去修改参数a，b

for times in range(1000):

sum1=0

sum2=0

for i in range(m):

sum1=sum1+diff(i)

sum2=sum2+diff(i)*X[i]

b=b-alpha*sum1

a=a-alpha*sum2

print ("参数b的值:%.2f"%b)

print ("参数a的值:%.2f"%a)

打印结果：
参数b的值:2.44
参数a的值:1.95

通过迭代方法求得的参数a，b与通过梯度求得的参数进行比较，发现他们非常接近，有异曲同工之妙！
这种方法也可看作只有一个神经元的神经网络，并且没有激活函数的这种。

16.7 迭代中用矩阵代替循环

这里利用梯度下降更新参数时采用一个循环，循环一般比较耗费资源，如果有一千万个数据，将需要循环一千万次，这是不可接受的，那么我们是否能不用循环？
当然可以，如果用采用矩阵方式，需要进行如下操作：
 把输入X,Y变为矩阵；
 把模型变为矩阵或向量；
 把式（1.13）变为矩阵与向量的点乘
 把式（1.14）变为矩阵的累加
具体代码实现如下：
（1）定义一个线性回归类，在这个类中初始化两个参数，并定义几个函数。

class LinearRegressGD(object):
    def __init__(self,eta=0.001,n_iter=4000):
        self.eta=eta
        self.n_iter=n_iter
    def fit(self,X,Y):
        #初始化参数
        self.w=np.zeros(1+X.shape[1])
        self.w[0]=0.1
        self.w[1]=0.2
        self.cost=[]
        #开始迭代循环
        for i in range(self.n_iter):
            output=self.net_input(X)
            errors=(Y-output)
            self.w[1] +=self.eta*X.T.dot(errors)
            self.w[0] +=self.eta*(errors).sum()
            cost1=(errors**2).sum()/2.0
            self.cost.append(cost1)
        return self
    def net_input(self,X):
        return np.dot(X,self.w[1])+self.w[0]
    def predict(self,X):
        return self.net_input(X)

class LinearRegressGD(object):

def __init__(self,eta=0.001,n_iter=4000):

self.eta=eta

self.n_iter=n_iter

def fit(self,X,Y):

#初始化参数

self.w=np.zeros(1+X.shape[1])

self.w[0]=0.1

self.w[1]=0.2

self.cost=[]

#开始迭代循环

for i in range(self.n_iter):

output=self.net_input(X)

errors=(Y-output)

self.w[1] +=self.eta*X.T.dot(errors)

self.w[0] +=self.eta*(errors).sum()

cost1=(errors**2).sum()/2.0

self.cost.append(cost1)

return self

def net_input(self,X):

return np.dot(X,self.w[1])+self.w[0]

def predict(self,X):

return self.net_input(X)

（2）把输入变为10x1矩阵，如何运行以上函数

import numpy as np
X=np.array([4,8,9,8,7,12,6,10,6,9]).reshape(10,1)
Y=np.array([9,20,22,15,17,23,18,25,10,20]).reshape(10,1)
lr=LinearRegressGD()
lr.fit(X,Y)
print("参数a的值:%.3f"%lr.w[1])
print("参数b的值:%.3f"%lr.w[0])

import numpy as np

X=np.array([4,8,9,8,7,12,6,10,6,9]).reshape(10,1)

Y=np.array([9,20,22,15,17,23,18,25,10,20]).reshape(10,1)

lr=LinearRegressGD()

lr.fit(X,Y)

print("参数a的值:%.3f"%lr.w[1])

print("参数b的值:%.3f"%lr.w[0])

打印结果：
参数a的值:1.995
参数b的值:2.130
（3）把迭代过程可视化

import matplotlib.pyplot as plt
plt.plot(range(1,lr.n_iter+1),lr.cost)
plt.ylabel('SSE')
plt.xlabel('Epoch')
plt.show()

import matplotlib.pyplot as plt

plt.plot(range(1,lr.n_iter+1),lr.cost)

plt.ylabel('SSE')

plt.xlabel('Epoch')

plt.show()

16.8 用Tensorflow架构实现自动求导,求参数

import numpy as np 
import tensorflow as tf  #导入tensorflow 
import matplotlib.pyplot as plt # import matplotlib
%matplotlib inline

##构造训练数据##
x_data=np.array([4,8,9,8,7,12,6,10,6,9])
y_data=np.array([9,20,22,15,17,23,18,25,10,20])
#可视化
plt.figure() # Create a new figure
plt.scatter(x_data,y_data) #Plot a scatter draw of the random datapoints
plt.plot (x_data, 2.0 + 2.0 * x_data) # Draw one line with the line function


##初始化权重变量#
w=tf.Variable(tf.random_uniform([1],-1.0,1.0)) 
b=tf.Variable(tf.zeros([1]))     
y=w*x_data+b


##计算损失函数
loss=tf.reduce_mean(tf.square(y-y_data))  #判断与正确值的差距
optimizer=tf.train.GradientDescentOptimizer(0.01) #根据差距进行反向传播修正参数
train=optimizer.minimize(loss) #建立训练器

init=tf.global_variables_initializer() #初始化TensorFlow训练结构
sess=tf.Session()  #建立TensorFlow训练会话
#sess = tf.InteractiveSession()  
sess.run(init)     #将训练结构装载到会话中
print('初始权重值w:',w.eval(session=sess))
for  step in range(1000): #循环训练100次
     sess.run(train)  #使用训练器根据训练结构进行训练
     if  step%100==0:  #每20次打印一次训练结果
        print(step,sess.run(w),sess.run(b)) #训练次数，W值，b值
        
print(sess.run(loss))        
print('最后权重值w:',w.eval(session=sess))
print('最后偏移量b:',b.eval(session=sess))
#可视化预测直线
plt.figure() # Create a new figure
plt.scatter(x_data,y_data)
plt.plot (x_data, sess.run(b) + x_data * sess.run(w))

import numpy as np

import tensorflow as tf #导入tensorflow

import matplotlib.pyplot as plt # import matplotlib

%matplotlib inline

##构造训练数据##

x_data=np.array([4,8,9,8,7,12,6,10,6,9])

y_data=np.array([9,20,22,15,17,23,18,25,10,20])

#可视化

plt.figure() # Create a new figure

plt.scatter(x_data,y_data) #Plot a scatter draw of the random datapoints

plt.plot (x_data, 2.0 + 2.0 * x_data) # Draw one line with the line function

##初始化权重变量#

w=tf.Variable(tf.random_uniform([1],-1.0,1.0))

b=tf.Variable(tf.zeros([1]))

y=w*x_data+b

##计算损失函数

loss=tf.reduce_mean(tf.square(y-y_data)) #判断与正确值的差距

optimizer=tf.train.GradientDescentOptimizer(0.01) #根据差距进行反向传播修正参数

train=optimizer.minimize(loss) #建立训练器

init=tf.global_variables_initializer() #初始化TensorFlow训练结构

sess=tf.Session() #建立TensorFlow训练会话

#sess = tf.InteractiveSession()

sess.run(init) #将训练结构装载到会话中

print('初始权重值w:',w.eval(session=sess))

for step in range(1000): #循环训练100次

sess.run(train) #使用训练器根据训练结构进行训练

if step%100==0: #每20次打印一次训练结果

print(step,sess.run(w),sess.run(b)) #训练次数，W值，b值

print(sess.run(loss))

print('最后权重值w:',w.eval(session=sess))

print('最后偏移量b:',b.eval(session=sess))

#可视化预测直线

plt.figure() # Create a new figure

plt.scatter(x_data,y_data)

plt.plot (x_data, sess.run(b) + x_data * sess.run(w))

运行结果
初始权重值w: [-0.99174666]
0 [ 3.35317731] [ 0.51469594]
100 [ 2.17306757] [ 0.62030703]
200 [ 2.14828968] [ 0.83054471]
300 [ 2.12670541] [ 1.01368761]
400 [ 2.10790277] [ 1.17322707]
500 [ 2.09152317] [ 1.31220567]
600 [ 2.07725477] [ 1.43327284]
700 [ 2.06482506] [ 1.53873765]
800 [ 2.05399752] [ 1.63060951]
900 [ 2.0445652] [ 1.71064162]
6.90384
最后权重值w: [ 2.03642535]
最后偏移量b: [ 1.77970874]

16.9 拓展

15.1 聚类概述

“物以类聚人以群分",现实生活中很多事物都存在类似现象，当然这些现象反映在数据上，就需要我们通过一定算法，找出这些类或族。聚类算法就是解决类似问题而提出的。
聚类就是按照某个特定标准(如距离准则)把一个数据集分割成不同的类或簇，使得同一个簇内的数据对象的相似性尽可能大，同时不在同一个簇中的数据对象的差异性也尽可能地大。即聚类后同一类的数据尽可能聚集到一起，不同数据尽量分离。
目前，有大量的聚类算法。而对于具体应用，聚类算法的选择取决于数据的类型、聚类的目的。如果聚类分析被用作描述或探查的工具，可以对同样的数据尝试多种算法，以发现数据可能揭示的结果。
主要的聚类算法可以划分为如下几类：划分方法、层次方法、基于密度的方法、基于网格的方法以及基于模型的方法
每一类中都存在着得到广泛应用的算法，例如：划分方法中的k-means聚类算法、层次方法中的层次聚类算法、基于模型方法中的高斯混合聚类算法等。
目前,聚类问题的研究不仅仅局限于上述的硬聚类，即每一个数据只能被归为一类，模糊聚类也是聚类分析中研究较为广泛的一个分支。模糊聚类通过隶属函数来确定每个数据隶属于各个簇的程度，而不是将一个数据对象硬性地归类到某一簇中。目前已有很多关于模糊聚类的算法被提出，如著名的高斯混合聚类等。
下图演示了K-Means进行聚类的迭代过程：

下图为高斯混合聚类迭代过程：

15.2 k-means模型

算法步骤：
（1）首先我们选择一些类/组，并随机初始化它们各自的中心点。
（2）计算每个数据点到中心点的距离，数据点距离哪个中心点最近就划分到哪一类中。
（3）重新计算每一类中心点作为新的中心点，各中心点求每个类中的平均值。
（4）重复以上步骤，直到每一类中心在每次迭代后变化不大为止。也可以多次随机初始化中心点，然后选择运行结果最好的一个。

优点：
计算简便
缺点：
我们必须提前知道数据有多少类/组。

15.3 简单实例

假定我们有如下8个点：
A1(2, 10) A2(2, 5) A3(8, 4) A4(5, 8) A5(7, 5) A6(6, 4) A7(1, 2) A8(4, 9)
现希望分成3个聚类（即k=3）
初始化选择 A1(2, 10), A4(5, 8) ，A7(1, 2)为聚类中心点，假设两点距离定义为ρ(a, b) = |x2 – x1| + |y2 – y1| . （当然也可以定义为其它格式，如欧氏距离）
第一步：选择3个聚类中，分别为A1，A4，A7

这些点的分布图如下：

图1
第二步：计算各点到3个类中心的距离，那个点里类中心最近，就把这个样本点
划归到这个类。选定3个类中心（即:A1，A4，A7），如下图：
图2
对A1点，计算其到每个cluster 的距离
A1->class1 = |2-2|+|10-10}=0
A1->class2 = |2-5|+|10-8|=5
A1->class3 = |2-1|+|10-2|=9
因此A1 属于cluster1，如下表：
(
按照类似方法，算出各点到聚类中心的距离，然后按照最近原则，把样本点放在那个族中。如下表：

根据距离最短原则，样本点的第一次样本划分，如下图：
图3
第三步：求出各类中样本点的均值，并以此为类的中心。
cluster1只有1个点，因此A1为中心点
cluster2的中心点为 ( (8+5+7+6+4)/5,(4+8+5+4+9)/5 )=（6,6）。注意：这个点并非样本点。
cluster3的中心点为( (2+1)/2, (5+2)/2 )= (1.5, 3.5)，
新族的具体请看下图中x点：

图4
第四步：计算各样本点到各聚类中心的距离，重复以上第二、第三步，把样本划分到新聚类中，如下图：
图5
持续迭代，直到前后两次迭代不发生变化为止，如下：
图6

15.4 简单实例用Python实现

（1）生成数据

from numpy import *
dataSet=array([[2,10],[2,5],[8,4],[5,8],[7,5],[6,4],[1,2],[4,9]])

1 2	from numpy import * dataSet=array([[2,10],[2,5],[8,4],[5,8],[7,5],[6,4],[1,2],[4,9]])

（2）创建距离函数

def distabs(vecA, vecB):
    '''
    # 计算两个向量的距离，用的是差的绝对值
    :param vecA:
    :param vecB: 
    :return: 距离
    '''
    return sum(abs(vecA - vecB))

def distabs(vecA, vecB):

'''

# 计算两个向量的距离，用的是差的绝对值

:param vecA:

:param vecB:

:return: 距离

'''

return sum(abs(vecA - vecB))

（3）手工选择3个聚类中心

def selectCent(dataSet, k):
    '''
    :param dataSet: 输入需聚类的数据集
    :param k: 划分的簇数
    :return: 返回 k 个簇的初始随机中心数组
    '''
    centroids = []
    for j in [0,3,6]:
        dataSet[j]
        centroids.append(dataSet[j])
    return array(centroids)

def selectCent(dataSet, k):

'''

:param dataSet: 输入需聚类的数据集

:param k: 划分的簇数

:return: 返回 k 个簇的初始随机中心数组

'''

centroids = []

for j in [0,3,6]:

dataSet[j]

centroids.append(dataSet[j])

return array(centroids)

（4）创聚类函数

def kMeans(dataSet, k,max_times=6, distMeas=distabs, createCent=selectCent):
    '''
    # K-mean 算法实现
    :param dataSet: 需聚类的数据集
    :param k: 要划分的簇数
    :param max_times:最高迭代次数 
    :param distMeas: 距离的度量算法，默认方法为 distabs（即差的绝对值距离）
    :param selectCent: 生成初始的质心
    :return: 返回 K 个簇的中心以及划分过后的数据集
    '''
    m = shape(dataSet)[0]

    # 创建一个数据划分空间
    clusterAssment = mat(zeros((m, 2)))
    # 获取初始簇中心点
    centroids = createCent(dataSet, k)
    clusterChanged = True
    while clusterChanged or max_times < 0:
        max_times -= 1
        clusterChanged = False
        for i in range(m):
            minDist = inf
            minIndex = -1
            for j in range(k):
                # 获取每个点到中心 k 的欧几里得距离
                distJI = distMeas(centroids[j, :], dataSet[i, :])
                # 选择最小距离
                if distJI < minDist:
                    minDist = distJI
                    minIndex = j
            # 判断中心点是否变化
            if clusterAssment[i, 0] != minIndex:
                clusterChanged = True
            clusterAssment[i, :] = minIndex, minDist ** 2

        # 重新计算中心
        for cent in range(k):
            # 获取这个簇中的所有点
            ptsInClust = dataSet[nonzero(clusterAssment[:, 0].A == cent)[0]]
            # 分配质心为平均值
            centroids[cent, :] = mean(ptsInClust, axis=0)
    return centroids, clusterAssment

def kMeans(dataSet, k,max_times=6, distMeas=distabs, createCent=selectCent):

'''

# K-mean 算法实现

:param dataSet: 需聚类的数据集

:param k: 要划分的簇数

:param max_times:最高迭代次数

:param distMeas: 距离的度量算法，默认方法为 distabs（即差的绝对值距离）

:param selectCent: 生成初始的质心

:return: 返回 K 个簇的中心以及划分过后的数据集

'''

m = shape(dataSet)[0]

# 创建一个数据划分空间

clusterAssment = mat(zeros((m, 2)))

# 获取初始簇中心点

centroids = createCent(dataSet, k)

clusterChanged = True

while clusterChanged or max_times < 0:

max_times -= 1

clusterChanged = False

for i in range(m):

minDist = inf

minIndex = -1

for j in range(k):

# 获取每个点到中心 k 的欧几里得距离

distJI = distMeas(centroids[j, :], dataSet[i, :])

# 选择最小距离

if distJI < minDist:

minDist = distJI

minIndex = j

# 判断中心点是否变化

if clusterAssment[i, 0] != minIndex:

clusterChanged = True

clusterAssment[i, :] = minIndex, minDist ** 2

# 重新计算中心

for cent in range(k):

# 获取这个簇中的所有点

ptsInClust = dataSet[nonzero(clusterAssment[:, 0].A == cent)[0]]

# 分配质心为平均值

centroids[cent, :] = mean(ptsInClust, axis=0)

return centroids, clusterAssment

（5）运行

kMeans(dataSet,3)

1	kMeans(dataSet,3)

运行结果
(array([[3, 9],
[7, 4],
[1, 3]]), matrix([[ 0., 4.],
[ 2., 9.],
[ 1., 1.],
[ 0., 9.],
[ 1., 1.],
[ 1., 1.],
[ 2., 1.],
[ 0., 1.]]))
这个运行结果与图6的结果一致。

15.5 简单实例用sklearn实现

（1）导入需要的库或模块

%matplotlib inline
import numpy as np      #科学计算包  
import matplotlib.pyplot as plt      #python画图包  
from matplotlib.font_manager import FontProperties
  
from sklearn.cluster import KMeans       #导入K-means算法包  

import matplotlib.font_manager as fm    ###便于中文显示
#myfont = fm.FontProperties(fname='/home/hadoop/anaconda3/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/simhei.ttf')
myfont = FontProperties(fname=r"c:\windows\fonts\simkai.ttf", size=14)

%matplotlib inline

import numpy as np #科学计算包

import matplotlib.pyplot as plt #python画图包

from matplotlib.font_manager import FontProperties

from sklearn.cluster import KMeans #导入K-means算法包

import matplotlib.font_manager as fm ###便于中文显示

#myfont = fm.FontProperties(fname='/home/hadoop/anaconda3/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/simhei.ttf')

myfont = FontProperties(fname=r"c:\windows\fonts\simkai.ttf", size=14)

（2）创建数据

X=np.array([[2,10],[2,5],[8,4],[5.0,8],[7,5],[6,4.0],[1,2],[4,9]])
y=np.array([0,1,1,1,1,2,2,2])

1 2	X=np.array([[2,10],[2,5],[8,4],[5.0,8],[7,5],[6,4.0],[1,2],[4,9]]) y=np.array([0,1,1,1,1,2,2,2])

（3）利用kmeans进行聚类，并把结果可视化

plt.figure(figsize=(12, 12))  
  
''''' 
centers:产生数据的聚类中心点，默认值3 
init:采用随机还是k-means++等（可使聚类中心尽可能的远的一种方法）
random_state:随机生成器的种子 

'''  
random_state = 10  
  
# 使用k-means聚类 

plt.subplot(221)  #在2图里添加子图1  
y_pred = KMeans(n_clusters=3, init='random',random_state=random_state).fit_predict(X) 
plt.scatter(X[:, 0], X[:, 1], c=y_pred)  
plt.title("使用kmeas聚类",fontproperties=myfont,size=12)   #加标题  
  
#使用kmeans++聚类
y_pred = KMeans(n_clusters=3, init='k-means++',random_state=random_state).fit_predict(X) 
  
plt.subplot(222)#在2图里添加子图2  
plt.scatter(X[:, 0], X[:, 1], c=y_pred)  
plt.title("使用kmeans++聚类",fontproperties=myfont,size=12)  

plt.show()

plt.figure(figsize=(12, 12))

'''''

centers:产生数据的聚类中心点，默认值3

init:采用随机还是k-means++等（可使聚类中心尽可能的远的一种方法）

random_state:随机生成器的种子

'''

random_state = 10

# 使用k-means聚类

plt.subplot(221) #在2图里添加子图1

y_pred = KMeans(n_clusters=3, init='random',random_state=random_state).fit_predict(X)

plt.scatter(X[:, 0], X[:, 1], c=y_pred)

plt.title("使用kmeas聚类",fontproperties=myfont,size=12) #加标题

#使用kmeans++聚类

y_pred = KMeans(n_clusters=3, init='k-means++',random_state=random_state).fit_predict(X)

plt.subplot(222)#在2图里添加子图2

plt.scatter(X[:, 0], X[:, 1], c=y_pred)

plt.title("使用kmeans++聚类",fontproperties=myfont,size=12)

plt.show()

15.6 简单实例用Tensorflow实现

这里需要用到很多tensorflow函数，大家可参考：
https://www.cnblogs.com/wuzhitj/p/6648563.html
(1)导入需要的库，初始化参数

import tensorflow as tf
import numpy as np
import time
import matplotlib
import matplotlib.pyplot as plt
from sklearn.datasets.samples_generator import make_blobs
#样本数
N=8
#族数
K=3
#最大迭代数
MAX_ITERS = 20
changed = True
iters = 0

import tensorflow as tf

import numpy as np

import time

import matplotlib

import matplotlib.pyplot as plt

from sklearn.datasets.samples_generator import make_blobs

#样本数

N=8

#族数

K=3

#最大迭代数

MAX_ITERS = 20

changed = True

iters = 0

（2）创建数据集，并可视化三个族中心

centers = [(2, 10.0), (5, 8), (1, 2)]
data=np.array([[2,10],[2,5],[8,4],[5.0,8],[7,5],[6,4.0],[1,2],[4,9]])
features=np.array([0,1,1,1,1,2,2,2])
fig, ax = plt.subplots()
#s-表示标记的大小
ax.scatter(np.asarray(centers).transpose()[0], np.asarray(centers).transpose()[1], marker = 'o', s = 250)
plt.show()

centers = [(2, 10.0), (5, 8), (1, 2)]

data=np.array([[2,10],[2,5],[8,4],[5.0,8],[7,5],[6,4.0],[1,2],[4,9]])

features=np.array([0,1,1,1,1,2,2,2])

fig, ax = plt.subplots()

#s-表示标记的大小

ax.scatter(np.asarray(centers).transpose()[0], np.asarray(centers).transpose()[1], marker = 'o', s = 250)

plt.show()

（3）可视化样本

fig, ax = plt.subplots()
ax.scatter(np.asarray(centers).transpose()[0], np.asarray(centers).transpose()[1], marker = 'o', s = 250)
ax.scatter(data.transpose()[0], data.transpose()[1], marker = 'o', s = 100, c = features, cmap=plt.cm.coolwarm )
plt.show()

fig, ax = plt.subplots()

ax.scatter(np.asarray(centers).transpose()[0], np.asarray(centers).transpose()[1], marker = 'o', s = 250)

ax.scatter(data.transpose()[0], data.transpose()[1], marker = 'o', s = 100, c = features, cmap=plt.cm.coolwarm )

plt.show()

（4）计算各样本的到各族中心距离

points=tf.Variable(data)
cluster_assignments = tf.Variable(tf.zeros([N], dtype=tf.int64))
centroids = tf.Variable(tf.slice(points.initialized_value(), [0,0], [K,2]))
sess = tf.Session()
sess.run(tf.global_variables_initializer())
sess.run(centroids)

#为计算每点对族中心的距离，使rep_centroids、rep_points都变为NxKx2矩阵
rep_centroids = tf.reshape(tf.tile(centroids, [N, 1]), [N, K, 2])
rep_points = tf.reshape(tf.tile(points, [1, K]), [N, K, 2])
sum_squares = tf.reduce_sum(tf.square(rep_points - rep_centroids),2)
#获取最小值的对应索引值
best_centroids = tf.argmin(sum_squares, 1)
#判断所有族中心是否不再变化
did_assignments_change = tf.reduce_any(tf.not_equal(best_centroids, cluster_assignments))

points=tf.Variable(data)

cluster_assignments = tf.Variable(tf.zeros([N], dtype=tf.int64))

centroids = tf.Variable(tf.slice(points.initialized_value(), [0,0], [K,2]))

sess = tf.Session()

sess.run(tf.global_variables_initializer())

sess.run(centroids)

#为计算每点对族中心的距离，使rep_centroids、rep_points都变为NxKx2矩阵

rep_centroids = tf.reshape(tf.tile(centroids, [N, 1]), [N, K, 2])

rep_points = tf.reshape(tf.tile(points, [1, K]), [N, K, 2])

sum_squares = tf.reduce_sum(tf.square(rep_points - rep_centroids),2)

#获取最小值的对应索引值

best_centroids = tf.argmin(sum_squares, 1)

#判断所有族中心是否不再变化

did_assignments_change = tf.reduce_any(tf.not_equal(best_centroids, cluster_assignments))

（5）定义函数，更新各族中心坐标

#定义函数，更新各族中心坐标
def bucket_mean(data, bucket_ids, num_buckets):
    total = tf.unsorted_segment_sum(data, bucket_ids, num_buckets)
    count = tf.unsorted_segment_sum(tf.ones_like(data), bucket_ids, num_buckets)
    return total / count
means = bucket_mean(points, best_centroids, K)
#确定执行依赖关系，先执行did_assignments_change，然后执行后续命令
with tf.control_dependencies([did_assignments_change]):
    do_updates = tf.group(centroids.assign(means),cluster_assignments.assign(best_centroids))

#定义函数，更新各族中心坐标

def bucket_mean(data, bucket_ids, num_buckets):

total = tf.unsorted_segment_sum(data, bucket_ids, num_buckets)

count = tf.unsorted_segment_sum(tf.ones_like(data), bucket_ids, num_buckets)

return total / count

means = bucket_mean(points, best_centroids, K)

#确定执行依赖关系，先执行did_assignments_change，然后执行后续命令

with tf.control_dependencies([did_assignments_change]):

do_updates = tf.group(centroids.assign(means),cluster_assignments.assign(best_centroids))

（6）可视化迭代过程

fig, ax = plt.subplots()
colourindexes=[2,1,3]
#循环停止条件是族中心不再变化而且循环次数不超过指定值
while changed and iters < MAX_ITERS:
    fig, ax = plt.subplots()
    iters += 1
    [changed, _] = sess.run([did_assignments_change, do_updates])
    [centers, assignments] = sess.run([centroids, cluster_assignments])
    ax.scatter(sess.run(points).transpose()[0], sess.run(points).transpose()[1], marker = 'o', s = 200, c = assignments, cmap=plt.cm.coolwarm )
    ax.scatter(centers[:,0],centers[:,1], marker = '^', s = 550, c = colourindexes, cmap=plt.cm.plasma)
    ax.set_title('Iteration ' + str(iters))
    plt.savefig("kmeans" + str(iters) +".png")
ax.scatter(sess.run(points).transpose()[0], sess.run(points).transpose()[1], marker = 'o', s = 200, c = assignments, cmap=plt.cm.coolwarm )
plt.show()

fig, ax = plt.subplots()

colourindexes=[2,1,3]

#循环停止条件是族中心不再变化而且循环次数不超过指定值

while changed and iters < MAX_ITERS:

fig, ax = plt.subplots()

iters += 1

[changed, _] = sess.run([did_assignments_change, do_updates])

[centers, assignments] = sess.run([centroids, cluster_assignments])

ax.scatter(sess.run(points).transpose()[0], sess.run(points).transpose()[1], marker = 'o', s = 200, c = assignments, cmap=plt.cm.coolwarm )

ax.scatter(centers[:,0],centers[:,1], marker = '^', s = 550, c = colourindexes, cmap=plt.cm.plasma)

ax.set_title('Iteration ' + str(iters))

plt.savefig("kmeans" + str(iters) +".png")

ax.scatter(sess.run(points).transpose()[0], sess.run(points).transpose()[1], marker = 'o', s = 200, c = assignments, cmap=plt.cm.coolwarm )

plt.show()

迭代1次就到达最佳结果，看来是要tensorflow效果不错！

15.7 改进

由于 K-means 算法的分类结果会受到初始点的选取而有所区别，因此有提出这种算法的改进: K-means++。其实这个算法也只是对初始点的选择有改进而已，其他步骤都一样。初始质心选取的基本思路就是，初始的聚类中心之间的相互距离要尽可能的远。整个算法的过程如下：
下面结合一个简单的例子说明K-means++是如何选取初始聚类中心的。数据集中共有8个样本，分布以及对应序号如下图所示：
图7
假设经过图7的步骤一后6号点被选择为第一个初始聚类中心，那在进行步骤二时每个样本的D(x)和被选择为第二个聚类中心的概率如下表所示：

其中的P(x)就是每个样本被选为下一个聚类中心的概率。最后一行的Sum是概率P(x)的累加和，用于轮盘法选择出第二个聚类中心。方法是随机产生出一个0~1之间的随机数，判断它属于哪个区间，那么该区间对应的序号就是被选择出来的第二个聚类中心了。例如1号点的区间为[0,0.2)，2号点的区间为[0.2, 0.525)。
从上表可以直观的看到第二个初始聚类中心是1号，2号，3号，4号中的一个的概率为0.9。而这4个点正好是离第一个初始聚类中心6号点较远的四个点。这也验证了K-means的改进思想：即离当前已有聚类中心较远的点有更大的概率被选为下一个聚类中心。可以看到，该例的K值取2是比较合适的。当K值大于2时，每个样本会有多个距离，需要取最小的那个距离作为D(x)。

Python技术交流与分享

分享技术平台

月度归档：2018年08月

16.1 目的

16.2 模型

16.3 问题

16.4 简单入手

16.5、用Python解方程求出a和b

16.6 用迭代方式求参数

16.7 迭代中用矩阵代替循环

16.8 用Tensorflow架构实现自动求导,求参数

16.9 拓展

15.1 聚类概述

15.2 k-means模型

15.3 简单实例

15.4 简单实例用Python实现

15.5 简单实例用sklearn实现

15.6 简单实例用Tensorflow实现

15.7 改进