主成分分析PCA+编程实现

主成分分析

  • 一个非监督的机器学习算法
  • 主要用于数据的降维
  • 通过降维,可以发现更便于人类理解的特征
  • 其他应用:可视化,去噪

对数据进行简化的一系列原因:

  • 使得数据集更易使用
  • 降低很多算法的计算开销
  • 去除噪声
  • 使得结果易懂

什么是主成分分析?

在PCA中,数据从原来的坐标系转换到了新的坐标系,新的坐标系的选择是由数据本身决定的。第一个新坐标轴选择的是原始数据中方差最大的方向,第二个新坐标轴的选择和第一个坐标轴正交且具有最大方差的方向。该过程一直重复,重复次数为原始数据中特征的数目,我们会发现,大部分方差都包含在最前面的几个新坐标轴中。因此,可以忽略余下的坐标轴,及对数据进行了降维处理。

PCA优缺点

  • 优点:降低数据的复杂性,识别最重要的多个特征
  • 缺点:不一定需要,且可能损失有用信息
  • 使用数据类型:数值型数据

举个例子理解一下,二维空间上有一些样本点如下,将二维降维至一维:

原始样本点(二维)映射到特征1上映射到特征2上

从图中可以看到,将样本点映射到特征1这种方案感觉会好一些,因为相较另一种方案,样本与样本之间的间距较大一些,也就是说点和点之间的区分度更大一些。那有没有更好地映射方案呢?来看看下面这中映射成一维的方案:

原始样本点(二维)找一根直线将所有的点都映射到这根直线上

在这种映射方案下,我们可以看到它将所有的点都映射到直线上,虽然线是斜的,我们也可以把它看成是一个轴,一个维度,可以看到所有的点在映射之后,可以发现它们与原来的样本点的分布并没有更大的差距,并且样本之间的距离(可区分度)也比无论映射到x轴还是y轴都要大。

那么我们需要思考的就是如何找到这个让样本间间距最大的轴?如何定义样本间的间距?通常情况下,我们使用方差(Variance),因为方差可以作为反映样本的疏密的一个指标。


这样上面的降维问题,就转化成为了找到一个轴,使得样本空间的所有点映射到这个轴后,方差最大的问题。

为了解决这个问题,在进行主成分分析之前,要做一些工作:

第一步:将样例的均值归为0(demean),所有样本减去整体的均值,可以看到处理之后,样本分布没有改变,只是坐标轴变了。

原始样本分布demean之后

demean之后样本均值为0,方差公式就可以做出一些改变,更加方便


第二步,求这个轴的方向w=(w1,w2)(两维的表达),使得我们所有的样本映射到w维度上方差最大,即


如何求图解

w为要映射的轴的方向,蓝色的点为原始的样本点,映射到w上,其实就是蓝色的那跟线段,其实也就是点乘的定义


由这个问题,扩展到N维,主成分分析的目标就是求w,使得映射到此维度后,方差最大


插曲,光看图我们会感觉到pca和线性回归有些相似,但需要注意他们之间的差别

PCA线性回归

横轴轴是特征

推导式子不同

虽然都是样本与一条直线的关系,但pca是距离时点垂直于直线

横轴特征,纵轴输出标记

推导式子不同

点垂直于特征

梯度上升法解决主成分分析问题


1.生成数据,进行demean

import numpy as np
import matplotlib.pyplot as plt
X=np.empty(shape=(100,2))
#从一个均匀分布[low,high)中随机采样,注意定义域是左闭右开,即包含low,不包含high.
X[:,0]=np.random.uniform(1,100,size=100)
X[:,1]=0.75*X[:,0]+3+np.random.normal(0,10,size=100)
plt.scatter(X[:,0],X[:,1])
plt.title('X')
plt.show()

def demean(X):
    #求每列的均值,也是每个特征维度上的均值
    return X-np.mean(X,axis=0)
X_demean=demean(X)
plt.scatter(X_demean[:,0],X_demean[:,1])
plt.title('demean(X)')
plt.show()

2.梯度上升法

def f(w,X):
    '''
    目标函数
    :param w:映射向量
    :param X: 样本数据
    :return:
    '''
    return np.sum(np.dot(X,w)**2)/len(X)

def df(w,X):
    '''
    目标函数求梯度
    :param w:
    :param X:
    :return:
    '''
    return 2/len(X)*np.dot(X.T,np.dot(X,w))

def direction(w):
    '''
    将向量w转化为单位向量
    :param w:
    :return:
    '''
    return w/np.linalg.norm(w)

def gradient_ascent(df,X,init_w,eta,max_iters=1e4,epsilon=1e-8):
    '''
    梯度上升
    :param df: 
    :param X: 
    :param init_w: 
    :param eta: 步长学习率
    :param max_iters: 最大迭代步数
    :param epsilon: 精度
    :return: 
    '''
    cur_iter=0
    w=direction(init_w)
    while cur_iter<max_iters:
        gradient=df(w,X)
        last_w=w
        w=w+eta*gradient
        w=direction(w)
        if(abs(f(w,X)-f(last_w,X))<epsilon):
            break
        cur_iter+=1
    return w

init_w=np.random.random(X.shape[1])
eta=0.01
w=gradient_ascent(df,X,init_w,eta)

plt.scatter(X_demean[:,0],X_demean[:,1])
plt.plot([-50,w[0]*80],[-50,w[1]*80],c='r')
plt.title('demean(X)')
plt.show()
求出第一主成分以后,如何求出下一个主成分?

将数据进行改变,将数据在第一个主成分上的分量去掉,即下图中的X',然后在新的数据上求第一主成分


# X2=np.empty(X.shape)
# for i in range(len(X)):
#     X2[i]=X[i]-np.dot(X[i],w)*w
X2=X-np.dot(X,w).reshape(-1,1)*w
plt.scatter(X2[:,0],X2[:,1])
plt.title('X2')
plt.show()
w2=gradient_ascent(df,X2,init_w,eta)
print(w2)
#垂直应该等于0
print(w.dot(w2))
将求前n个主成分,封装成函数
def first_n_components(n,X,eta=0.01,n_iter=1e4,epsilon=1e-8):
    x_pca=X.copy()
    x_pca=demean(x_pca)
    res=[]
    for i in range(n):
        init_w=np.random.random(x_pca.shape[1])
        w=gradient_ascent(df,x_pca,init_w,eta)
        res.append(w)
        x_pca=x_pca-np.dot(x_pca,w).reshape(-1,1)*w
    return res

print(first_n_components(2,X))
[array([ 0.75823566,  0.65198059]), array([ 0.65198379, -0.7582329 ])]

高维数据向低维数据映射

样本数据(m行n列)X的前k个主成分(k行,n列)

高维数据向低维数据映射,其实就是求


反之,从低维映射到高维也是可以的


import numpy as np
class PCA:
    def __init__(self,n_components):
        '''
        :param n_components: 指定几个主成分
        '''
        self.n_components=n_components
        self.components_=None

    def __repr__(self):
        return 'pca%d个主成分'%self.n_components

    def fit(self,X,eta=0.01,n_iter=1e4):
        '''根据训练数据集X获得数据的均值和方差'''
        def demean(X):
            # 求每列的均值,也是每个特征维度上的均值
            return X - np.mean(X, axis=0)

        def f(w, X):
            '''
            目标函数
            :param w:映射向量
            :param X: 样本数据
            :return:
            '''
            return np.sum(X.dot(w) ** 2) / len(X)

        def df(w, X):
            '''
            目标函数求梯度
            :param w:
            :param X:
            :return:
            '''
            return 2 / len(X) * np.dot(X.T, np.dot(X, w))

        def direction(w):
            '''
            将向量w转化为单位向量
            :param w:
            :return:
            '''
            return w / np.linalg.norm(w)

        def gradient_ascent(df, X, init_w, eta=eta, max_iters=n_iter, epsilon=1e-8):
            '''
            梯度上升
            :param df:
            :param X:
            :param init_w:
            :param eta: 步长学习率
            :param max_iters: 最大迭代步数
            :param epsilon: 精度
            :return:
            '''
            w = direction(init_w)
            cur_iter = 0

            while cur_iter < max_iters:
                gradient = df(w, X)
                last_w = w
                w = w + eta * gradient
                w = direction(w)
                if (abs(f(w, X) - f(last_w, X)) < epsilon):
                    break
                cur_iter += 1
            return w

        x_pca=demean(X)
        self.components_=np.empty(shape=(self.n_components,X.shape[1]))
        for i in range(self.n_components):
            init_w=np.random.random(x_pca.shape[1])
            w=gradient_ascent(df,x_pca,init_w)
            self.components_[i,:]=w
            x_pca=x_pca-np.dot(x_pca,w).reshape(-1,1)*w
        return self

    def transform(self,X):
        '''将给定的X,映射到各个主成分分量中'''
        return np.dot(X,np.transpose(self.components_))

    def inverse_transform(self,X):
        '''将给定的X,反向映射回原来的特征空间'''
        return np.dot(X,self.components_)

import matplotlib.pyplot as plt
X=np.empty(shape=(100,2))
#从一个均匀分布[low,high)中随机采样,注意定义域是左闭右开,即包含low,不包含high.
X[:,0]=np.random.uniform(1,100,size=100)
X[:,1]=0.75*X[:,0]+3+np.random.normal(0,10,size=100)
pca=PCA(n_components=2)
pca.fit(X)
print(np.transpose(pca.components_))
#降维
pca=PCA(n_components=1)
pca.fit(X)
x_reduction=pca.transform(X)
print(x_reduction.shape)

#恢复原始维度
x_restore=pca.inverse_transform(x_reduction)
print(x_restore.shape)
plt.scatter(X[:,0],X[:,1],c='b')
plt.scatter(x_restore[:,0],x_restore[:,1],c='g')
plt.show()

SKlearn中的PCA

import numpy as np
import matplotlib.pyplot as plt
X=np.empty(shape=(100,2))
#从一个均匀分布[low,high)中随机采样,注意定义域是左闭右开,即包含low,不包含high.
X[:,0]=np.random.uniform(1,100,size=100)
X[:,1]=0.75*X[:,0]+3+np.random.normal(0,10,size=100)

from sklearn.decomposition import PCA
pca=PCA(n_components=1)
pca.fit(X)
print(pca)
print(pca.components_)
x_reduction=pca.transform(X)
print(x_reduction.shape)
x_restore=pca.inverse_transform(x_reduction)
plt.scatter(X[:,0],X[:,1],c='b')
plt.scatter(x_restore[:,0],x_restore[:,1],c='g')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
digits=datasets.load_digits()
x=digits.data
y=digits.target

from sklearn.model_selection import train_test_split
x_train,x_test,y_train,y_test=train_test_split(x,y,random_state=666)
print(x_train.shape)

from sklearn.neighbors import KNeighborsClassifier
knn_clf=KNeighborsClassifier()
knn_clf.fit(x_train,y_train)
print(knn_clf.score(x_test,y_test))

from sklearn.decomposition import PCA
pca=PCA(n_components=2)
pca.fit(x_train)
x_train_reduction=pca.transform(x_train)
x_test_reduction=pca.transform(x_test)

knn_clf=KNeighborsClassifier()
knn_clf.fit(x_train_reduction,y_train)
print(knn_clf.score(x_test_reduction,y_test))
#解释各轴方差相应的比例
print(pca.explained_variance_ratio_)

pca=PCA(n_components=64)
pca.fit(x_train)
print(pca.explained_variance_ratio_)
plt.plot([i for i in range(64)],[np.sum(pca.explained_variance_ratio_[:i+1]) for i in range(64)])
plt.show()

pca=PCA(0.95)
pca.fit(x_train)
print(pca.n_components_)
x_train_reduction=pca.transform(x_train)
x_test_reduction=pca.transform(x_test)
knn_clf=KNeighborsClassifier()
knn_clf.fit(x_train_reduction,y_train)
print(knn_clf.score(x_test_reduction,y_test))

#降到2维进行可视化
pca=PCA(n_components=2)
pca.fit(x)
x_reduction=pca.transform(x)
for i in range(10):
    plt.scatter(x_reduction[y==i,0],x_reduction[y==i,1],alpha=0.5)
plt.show()

PCA+MNIST

import numpy as np
from sklearn.datasets import fetch_mldata
mnist=fetch_mldata("MNIST original",data_home='./mnist')
print(mnist)
x=mnist['data']
y=mnist['target']
x_train=np.array(x[:60000],dtype=float)
y_train=np.array(y[:60000],dtype=float)
x_test=np.array(x[60000:],dtype=float)
y_test=np.array(y[60000:],dtype=float)

from sklearn.neighbors import KNeighborsClassifier
knn_clf=KNeighborsClassifier()
knn_clf.fit(x_train,y_train)
print(knn_clf.score(x_test,y_test))


#pca降维
from sklearn.decomposition import PCA
pca=PCA(0.9)
pca.fit(x_train)
x_train_reduction=pca.transform(x_train)
knn_clf.fit(x_train_reduction,y_train)
x_test_reduction=pca.transform(x_test)
print(knn_clf.score(x_test_reduction,y_test))
{'COL_NAMES': ['label', 'data'], 'data': array([[0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       ..., 
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0]], dtype=uint8), 'target': array([ 0.,  0.,  0., ...,  9.,  9.,  9.]), 'DESCR': 'mldata.org dataset: mnist-original'}
0.9688
0.9728

可以看到在降维之后,不仅节省了计算开销,正确率也上升了,一箭双雕。


使用PCA对数据进行降噪

import numpy as np
from sklearn import datasets
import matplotlib.pyplot as plt
digits=datasets.load_digits()
x=digits.data
y=digits.target
#生成噪声数据
noise_digits=x+np.random.normal(0,4,size=x.shape)
#y=0的数字10个
example_digits=noise_digits[y==0,:][:10]
print(example_digits.shape)
#将另外的1-9标签各类图片10张,压入example_digits中
for num in range(1,10):
    x_num=noise_digits[y==num,:][:10]
    example_digits=np.vstack([example_digits,x_num])
print(example_digits.shape)
#画图
def plot_digits(data,title):
    fig,axes=plt.subplots(10,10,figsize=(10,10),
                         subplot_kw={'xticks':[],'yticks':[]},gridspec_kw=dict(hspace=0.1,wspace=0.1))
    for i,ax in enumerate(axes.flat):
        ax.imshow(data[i].reshape(8,8),
                  cmap='binary',
                  interpolation='nearest',
                  clim=(0,16))
    plt.title(title)
    plt.show()
plot_digits(example_digits,'降噪前')
#因为噪声较多,只选取50%的信息
from sklearn.decomposition import PCA
pca=PCA(0.5)
pca.fit(noise_digits)
print(pca.n_components_)
components=pca.transform(example_digits)
filtered_digits=pca.inverse_transform(components)
plot_digits(filtered_digits,'降噪后')
降噪前降噪后
特征脸

import numpy as np
from sklearn.datasets import fetch_lfw_people
import matplotlib.pyplot as plt
faces=fetch_lfw_people()


#对所有的样本进行一个随机的排列
random_indexes=np.random.permutation(len(faces.data))
x=faces.data[random_indexes]

example_faces=x[:36,:]
print(example_faces.shape)
def plot_faces(data):
    fig,axes=plt.subplots(6,6,figsize=(10,10),
                         subplot_kw={'xticks':[],'yticks':[]},gridspec_kw=dict(hspace=0.1,wspace=0.1))
    for i,ax in enumerate(axes.flat):
        #图片w*h=62,47
        ax.imshow(data[i].reshape(62,47))
    plt.show()
plot_faces(example_faces)

from sklearn.decomposition import PCA
pca=PCA(svd_solver='randomized')
pca.fit(x)
plot_faces(pca.components_[:36,:])
降维前降维后











评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值