plink 进行PCA分析

Posted

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了plink 进行PCA分析相关的知识,希望对你有一定的参考价值。

参考技术A 本次使用train.vcf.gz 作为例子

得到3个以.map, .nosex, .ped结尾的文件。

得到2个以.bim, .bed结尾的文件。

得到2个以.eigenval, .eigenvec结尾的文件;其中.eigenval代表每个PCA所占的比重,另外一个记录特征向量,用于坐标轴的绘制

结果

** 若想分析部分样本,则可以使用--remove参数,后接一个文件,其格式为: 第一列:群体编号, 第二列:样本名称,在这个例子中

对图像进行主成分分析(PCV.tools.pca.pca)

1 引言

  1.1 维度灾难

    分类为例:如最近邻分类方法(基本思想:以最近的格子投票分类)

    问题:当数据维度增大,分类空间爆炸增长。如图1所示,

      技术分享图片

 

                      图1 维度增加示意图

 

  1.2 解决方法

    缓解维度遭难的一个重用途径是降维。降维是通过某种数学变换,将原始高维属性空间转换为一个低维“子空间”,在这个子空间中样本

  密度大幅度提高,距离计算也变得更为容易。

 

  1.3 降维的可行性

    数据样本虽然是高维的,但与我们关心的也许仅是某个低维分布,因而可以对数据进行有效的降维。

 

  1.4 主成分分析

  1.4.1 概念

    主成分分析(Principal Component Analysys,简称PCV),是将具有相关性的高维指标,通过线性变换,化为少数相互独立的低维综合指标的一种

  多元统计方法。 又称为主分量分析。使用PCV降维时应该满足:

               (1)使降维后样本的方差尽可能大。

               (2)使降维后数据的均方误差尽可能小。

  1.4.2 基本计算步骤

    (1)计算给定样本{xn},n=1,2,...,N的均值mean_x和方差S;

    (2)计算S的特征向量与特征值,X = UΛUT

    (3)将特征值从小到大排列,前M个特征值λ1,λ1,...,λM所对应的特征向量u1,u2,...,uM构成投影矩阵。

 

2 矩阵特征值及特征向量求解

  2.1 定义

    A为n阶矩阵,若数λ和n维非0列向量x满足Ax=λx,那么数λ称为A的特征值,x称为A的对应于特征值λ的特征向量。

    式Ax=λx也可写成( A-λE)x=0,并且|λE-A|叫做A 的特征多项式。

    |λE-A| = 0,称为A的特征方程(一个齐次线性方程组),求解特征值的过程其实就是求解特征方程的解。

 

  2.2 特征值分解

    特征值分解是将一个矩阵分解成下面的形式:

             技术分享图片

    其中Q是这个矩阵A的特征向量组成的矩阵Σ?=?diag(λ1,?λ2,?...,?λn)是一个对角阵,每一个对角线上的元素就是一个特征值。

 

3 奇异值分解(SVD)

  特征值分解只适用于方阵。而在现实的世界中,我们看到的大部分矩阵都不是方阵,我们怎样才能像描述特征值一样描述这样一般矩阵的重要特征呢?奇异值分解(Singular Value Decomposition, 简称SVD)可用来解决这个问题。SVD比特征值分解的使用范围更广,是一个适用于任一矩阵分解的方法。

  技术分享图片

 

   假设X是一个m * n的矩阵,那么得到的U是一个m * m的方阵(里面的向量是正交的,U里面的向量称为左奇异向量),

                   Σ是一个m * n的实数对角矩阵(对角线以外的元素都是0,对角线上的元素称为奇异值),

                   VT(V的转置)是一个n * n的矩阵,里面的向量也是正交的,V里面的向量称为右奇异向量,

  从图片来反映几个相乘的矩阵的大小可得下面的图片,

         技术分享图片

  将一个矩阵X的转置 XT * X,将会得到 XTX 是一个方阵,我们用这个方阵求特征值可以得到:

              (ATA)vi = λivi

     这里得到的V,就是我们上面的右奇异向量。这里的σi 就是就是上面说的奇异值λi,ui就是上面说的左奇异向量。常见的做法是将奇异值由大而小排列。如此Σ便能由M唯一确定了。

  奇异值σ跟特征值类似,在矩阵Λ中也是从大到小排列,而且σ的减少特别的快,在很多情况下,前10%甚至1%的奇异值的和就占了全部的奇异值之和的99%以上了。也就是说,我们也可以用前r大的奇异值来近似描述矩阵,这里定义一下部分奇异值分解

                   技术分享图片

    r是一个远小于m、n的数,这样矩阵的乘法看起来像是下面的样子:

          技术分享图片

   右边的三个矩阵相乘的结果将会是一个接近于X的矩阵,在这儿,r越接近于n,则相乘的结果越接近于A。而这三个矩阵的面积之和(在存储观点来说,矩阵面积越小,存储量就越小)要远远小于原始的矩阵A,我们如果想要压缩空间来表示原矩阵A,我们存下这里的三个矩阵:U、Λ、V就好了。

 

4 特征分解与奇异值分解

  4.1 低秩近似(Low-rank Approximation)

    技术分享图片

  

  4.2 关联

    技术分享图片

 

5 pca(matrix)函数代码

 1 def pca(X):
 2     """    主成分分析(PCA)
 3         输入: X, (每行存储训练数据(一张图像)的扁平数组)
 4         返回: 投影矩阵(V,重要维度优先)、方差(S)和均值(mean_X)
 5     """
 6     
 7     # 获得维度:num_data,dim = 图片数量 *  图片大小       ------>图片大小 = 长*宽 
 8     num_data,dim = X.shape
 9     
10     # 中心数据
11     mean_X = X.mean(axis=0)    # 每列求平均值 得到 1 * dim
12     X = X - mean_X             # X.shape = num_data * dim
13     
14     if dim > num_data:
15         # PCA - 紧凑的使用技巧
16         M = dot(X,X.T)      # 协方差矩阵
17         e,EV = linalg.eigh(M)  # e,EV = M的特征值,特征向量
18         tmp = dot(X.T,EV).T   # 使用紧凑技巧?
19         V = tmp[::-1]       # 为了得到最后的特征向量
20         S = sqrt(e)[::-1]     # 反向,因为特征值是递增的
21         for i in range(V.shape[1]):   # V.shape[1] = ???
22             V[:,i] /= S
23     else:
24         # PCA - 使用奇异值分解(SVD)
25         U,S,V = linalg.svd(X)
26         V = V[:num_data]      # 只有返回前num_data个才有意义       但是我认为此处不应该为V[:num_data]而应该为 V[:dim]  ------>因为num_data >= dim
27     
28     # 返回n the projection matrix, the variance and the mean
29     return V,S,mean_X

 

6 简单实例说明PCV实现的过程

6.1 列数大于行数时

 1 # 初始化矩阵X  (列数大于行数时)
 2 X =  [[1 0 1 0 0 0]
 3         [0 1 0 0 0 0]
 4         [1 1 0 0 0 0]
 5         [1 0 0 1 1 0]
 6         [0 0 0 1 0 1]]
 7 #计算每列的平均值
 8 mean_X =X.mean(axis=0)= [0.6 0.4 0.2 0.4 0.2 0.2]
 9 # 计算中心数据
10 X =X - mean_X =  [[ 0.4 -0.4  0.8 -0.4 -0.2 -0.2]
11  [-0.6  0.6 -0.2 -0.4 -0.2 -0.2]
12  [ 0.4  0.6 -0.2 -0.4 -0.2 -0.2]
13  [ 0.4 -0.4 -0.2  0.6  0.8 -0.2]
14  [-0.6 -0.4 -0.2  0.6 -0.2  0.8]]
15 # 计算协方差矩阵
16 M =dot(X, X.T)=[[ 1.20000000e+00 -4.00000000e-01  5.55111512e-17 -2.00000000e-01
17   -6.00000000e-01]
18  [-4.00000000e-01  1.00000000e+00  4.00000000e-01 -8.00000000e-01
19   -2.00000000e-01]
20  [ 5.55111512e-17  4.00000000e-01  8.00000000e-01 -4.00000000e-01
21   -8.00000000e-01]
22  [-2.00000000e-01 -8.00000000e-01 -4.00000000e-01  1.40000000e+00
23    0.00000000e+00]
24  [-6.00000000e-01 -2.00000000e-01 -8.00000000e-01  0.00000000e+00
25    1.60000000e+00]]
26 # 计算协方差的特征值和特征向量
27 e, EV = numpy.linalg.eigh(M) 
28 e =  [-6.04093951e-16  2.66097365e-01  1.17478886e+00  2.00000000e+00
29   2.55911378e+00]
30 EV =  [[ 4.47213595e-01 -1.38203855e-01  6.92423635e-01  5.00000000e-01
31   -2.26824169e-01]
32  [ 4.47213595e-01 -6.14411448e-01 -1.73966345e-01 -5.00000000e-01
33   -3.77139608e-01]
34  [ 4.47213595e-01  6.98980014e-01 -3.02870626e-01  4.38587603e-16
35   -4.68717745e-01]
36  [ 4.47213595e-01 -2.11286151e-01 -5.40988322e-01  5.00000000e-01
37    4.61183041e-01]
38  [ 4.47213595e-01  2.64921441e-01  3.25401658e-01 -5.00000000e-01
39    6.11498480e-01]]
40 #行之间反转
41 tmp = dot(X.T, EV).T =[[-3.33066907e-16 -3.33066907e-16  4.16333634e-16 -3.33066907e-16
42   -3.88578059e-16  1.66533454e-16]
43  [ 3.49490007e-01  8.45685660e-02 -1.38203855e-01  5.36352895e-02
44   -2.11286151e-01  2.64921441e-01]
45  [-1.51435313e-01 -4.76836971e-01  6.92423635e-01 -2.15586664e-01
46   -5.40988322e-01  3.25401658e-01]
47  [ 1.00000000e+00 -5.00000000e-01  5.00000000e-01 -7.21644966e-16
48    5.00000000e-01 -5.00000000e-01]
49  [-2.34358872e-01 -8.45857353e-01 -2.26824169e-01  1.07268152e+00
50    4.61183041e-01  6.11498480e-01]]
51 #默认返回的是按特征值升序,上下调转一下
52 V =tmp[::-1] = [[-2.34358872e-01 -8.45857353e-01 -2.26824169e-01  1.07268152e+00
53    4.61183041e-01  6.11498480e-01]
54  [ 1.00000000e+00 -5.00000000e-01  5.00000000e-01 -7.21644966e-16
55    5.00000000e-01 -5.00000000e-01]
56  [-1.51435313e-01 -4.76836971e-01  6.92423635e-01 -2.15586664e-01
57   -5.40988322e-01  3.25401658e-01]
58  [ 3.49490007e-01  8.45685660e-02 -1.38203855e-01  5.36352895e-02
59   -2.11286151e-01  2.64921441e-01]
60  [-3.33066907e-16 -3.33066907e-16  4.16333634e-16 -3.33066907e-16
61   -3.88578059e-16  1.66533454e-16]]
62 #sigma = 特征值开根号
63 #使特征值降序排序
64 S =sqrt(e)[::-1]= [1.59972303 1.41421356 1.08387677 0.51584626        nan]
65 #特征向量归一化
66 for i in range(V.shape[1]):
67      V[:, i] /= S
68 -----------------------------------------------------------------------------     
69 V =  [[-1.46499655e-01 -5.28752375e-01 -1.41789650e-01  6.70542025e-01
70    2.88289305e-01  3.82252720e-01]
71  [ 7.07106781e-01 -3.53553391e-01  3.53553391e-01 -5.10280049e-16
72    3.53553391e-01 -3.53553391e-01]
73  [-1.39716356e-01 -4.39936516e-01  6.38839814e-01 -1.98903298e-01
74   -4.99123458e-01  3.00220160e-01]
75  [ 6.77508074e-01  1.63941415e-01 -2.67916753e-01  1.03975338e-01
76   -4.09591321e-01  5.13566659e-01]
77  [            nan             nan             nan             nan
78               nan             nan]]
79 S =  [1.59972303 1.41421356 1.08387677 0.51584626        nan]

 

6.2 列数不大于行数时

 1  1 #初始化数组X
 2  2 X =  [[ 1  3]
 3  3  [ 4  7]
 4  4  [ 3  2]
 5  5  [ 1 23]]
 6  6 num_data, dim = X.shape
 7  7 X.shape= (4, 2)
 8  8 U, S, V = linalg.svd(X)
 9  9 U.shape, S.shape, V.shape=(4, 4) (2,) (2, 2)
10 10 V = [[ 0.10462808  0.99451142]
11 11  [ 0.99451142 -0.10462808]]
12 12 V = V[:dim]                                        #V = V[:dim]  or V =V[:num_data]    ???   注意:此时num_data > dim,所以我认为应该是V[:dim]
13 13 V[dim-1] =V[ 1 ]= [ 0.99451142 -0.10462808]
14 14 #计算结果
15 15 U = [[ 0.12635702  0.149642   -0.4064109  -0.89245244]
16 16  [ 0.30196809  0.71358512 -0.49810453  0.38923441]
17 17  [ 0.09422707  0.60994996  0.75324035 -0.22740114]
18 18  [ 0.94019702 -0.31042648  0.13910799  0.0177182 ]]
19 19 S =  [24.43997403  4.54836997]
20 20 V =  [[ 0.10462808  0.99451142]
21 21  [ 0.99451142 -0.10462808]]

 

7 环境配置

 7.1 需要准备的安装包

  下载python(x,y)2.7.x: http://www.softpedia.com/get/Programming/Other-Programming-Files/Python-x-y.shtml

  下载PVC包:https://github.com/willard-yuan/pcv-book-code

 

7.2 安装Python(x, y)

  在Windows下,推荐安装Python(x,y)2.7.x。Python(x,y)2.7.x是一个库安装包,除了Python自身外,还包含了第三方库,下面是安装Python(x, y)时的界面:

 

                                                        技术分享图片

    选择Full模式,自动安装的包见下图,

                                    技术分享图片

从上面第二幅图可以看出,pythonxy不仅包含了SciPy、NumPy、PyLab、OpenCV、MatplotLib,还包含了机器学习库scikits-learn。 为避免出现运行实例时出现的依赖问题,译者建议将上面的库全部选上,也就是选择“full”(译者也是用的全部安装的方式进行后面的实验的)。安装完成后,为验证安装是否正确,可以在Python shell里确认一下OpenCV是否已安装来进行验证,在Python Shell里输入下面命令:

1 from cv2 import __version__
2 __version__

 

自动安装的pyopencv有问题。需要你把opencv卸载后重新安装就可以了。 重新安装完成后,再次输入上面命令,就可以看到OpenCV的版本信息,则说明python(x,y)已安装正确。

  

  7.2 安装PCV库

  7.2.1 安装twine

    cd 到对应python的scripts目录中,执行以下命令

1 pip install twine  
或者
2 python3 -m pip install twine

  7.2.2 PVC的安装

    PCV库是一个第三方库。安装PVC之前必须安装twine,否则安装会出错。假设你已从上面网址上下载了中译版源码,从Windows cmd终端进入PCV所在目录:

1 cd PCV
2 python setup.py install

 

    运行上面命令,即可完成PCV库的安装。为了验证PCV库是否安装成功,在运行上面命令后,可以打开Python自带的Shell,在Shell输入:

 

1 import PCV

 

    如果未报错,则表明你已成功安装了该PCV库。

8 图像处理案案例

8.1 图像处理测试代码

 1 # -*- coding: utf-8 -*-
 2 import pickle
 3 from PIL import Image
 4 from numpy import *
 5 from pylab import *
 6 from PCV.tools import imtools, pca
 7 
 8 # 获取图像列表和他们的尺寸
 9 imlist = imtools.get_imlist(../data/fontimages/a_thumbs)  # fontimages.zip is part of the book data set
10 im = array(Image.open(imlist[0]))  # open one image to get the size
11 m, n = im.shape[:2]  # get the size of the images
12 imnbr = len(imlist)  # get the number of images
13 print "The number of images is %d" % imnbr
14 
15 # Create matrix to store all flattened images
16 immatrix = array([array(Image.open(imname)).flatten() for imname in imlist], f)
17 
18 # PCA降维
19 V, S, immean = pca.pca(immatrix)
20 
21 # 保存均值和主成分
22 f = open(../data/fontimages/font_pca_modes.pkl, wb)
23 pickle.dump(immean, f)
24 pickle.dump(V, f)
25 f.close()
26 
27 # Show the images (mean and 7 first modes)
28 # This gives figure 1-8 (p15) in the book.
29 figure()
30 gray()
31 subplot(2, 4, 1)
32 axis(off)
33 imshow(immean.reshape(m, n))
34 for i in range(7):
35     subplot(2, 4, i+2)
36     imshow(V[i].reshape(m, n))
37     axis(off)
38 show()

 

8.2 文件目录相对位置如下

data与执行文件的相对位置:

         技术分享图片

 

待处理图片.jpg的位置

      技术分享图片

8.3 执行结果

    技术分享图片

参考资料:http://yongyuan.name/pcvwithpython/installation.html

 




以上是关于plink 进行PCA分析的主要内容,如果未能解决你的问题,请参考以下文章

案例scatterplot3d 绘制 PCA_3D_plot

gwas分析时数据格式转换

使用 PCA 进行图像分析/特征提取

主成分分析(PCA)

使用 Python 进行 PCA 分析和绘图

PCA图像数据降维及重构误差分析实战并使用TSNE进行异常数据可视化分析