[python]用python实现的pca算法

pca算法用于原始数据维数较高时对数据进行降维

关于pca算法的学习，有一篇分析特别详细的论文http://www.cs.otago.ac.nz/cosc453/student_tutorials/principal_components.pdf  

比较好的中文总结：http://www.cnblogs.com/jerrylead/archive/2011/04/18/2020209.html

在使用python编写pca时需要使用python 的numpy，用这个算法库计算矩阵运算非常方便

由于没有接触过numpy以及网上介绍较少，代码是参考同学的代码写的。

#########################################################################
# File Name: npy.py
# Author: bbezxcy
# mail: [e-mail]
# Created Time: Thu 10 Jul 2014 06:45:16 PM CST
#########################################################################
#coding=utf-8
from numpy import *
def pca(mat, lenth):
	meanval = mean(mat, axis = 0)
	rmmeanMat = mat - meanval
	covMat = cov(rmmeanMat,rowvar = 0)
	eigval,eigvec =  linalg.eig(mat(covMat))
	tfMat =eigvec[0:lenth,:]
	finalData = rmmeanMat*tfMat 
	recoMat = finalData * tfMat.T + meanval
	return finalData,recoMat

 

 

编写代码过程中，为了方便，将自己的vimrc参考学弟的（http://noclyt.com/blog/?p=88)加上了一些功能（例如f5一键编译以及C艹用f8自动调用gdb等等），在这里给分享出来。

"编码设置
set enc=utf-8
set fencs=utf-8,ucs-bom,shift-jis,gb18030,gbk,gb2312,cp936
"语言设置
set langmenu=zh_CN.UTF-8
set helplang=cn
if has("syntax")  
  syntax on  
endif  
""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""  
" 显示相关    
""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""  
"set shortmess=atI   " 启动的时候不显示那个援助乌干达儿童的提示    
"winpos 5 5          " 设定窗口位置    
"set lines=40 columns=155    " 设定窗口大小    
"set nu              " 显示行号    
set go=             " 不要图形按钮    
"color asmanian2     " 设置背景主题    
set guifont=Courier_New:h10:cANSI   " 设置字体    
"syntax on           " 语法高亮    
autocmd InsertLeave * se nocul  " 用浅色高亮当前行    
autocmd InsertEnter * se cul    " 用浅色高亮当前行    
"set ruler           " 显示标尺    
set showcmd         " 输入的命令显示出来，看的清楚些    
"set cmdheight=1     " 命令行（在状态行下）的高度，设置为1    
"set whichwrap+=<,>,h,l   " 允许backspace和光标键跨越行边界(不建议)    
"set scrolloff=3     " 光标移动到buffer的顶部和底部时保持3行距离    
set novisualbell    " 不要闪烁(不明白)    
set statusline=%F%m%r%h%w\ [FORMAT=%{&ff}]\ [TYPE=%Y]\ [POS=%l,%v][%p%%]\ %{strftime(\"%d/%m/%y\ -\ %H:%M\")}   "状态行显示的内容    
set laststatus=1    " 启动显示状态行(1),总是显示状态行(2)    
set foldmethod=manual   " 手动折叠    
"set background=dark "背景使用黑色   
set nocompatible  "去掉讨厌的有关vi一致性模式，避免以前版本的一些bug和局限    
"colorscheme desert  
" 显示中文帮助  
if version >= 603  
    set helplang=cn  
    set encoding=utf-8  
endif  
" 设置配色方案  
"colorscheme murphy  
"字体   
"if (has("gui_running"))   
"   set guifont=Bitstream\ Vera\ Sans\ Mono\ 10   
"endif   
set fencs=utf-8,ucs-bom,shift-jis,gb18030,gbk,gb2312,cp936  
set termencoding=utf-8  
set encoding=utf-8  
set fileencodings=ucs-bom,utf-8,cp936  
set fileencoding=utf-8  
autocmd BufNewFile *.py,*.cpp,*.[ch],*.sh,*.java exec ":call SetTitle()"   
""定义函数SetTitle，自动插入文件头   
func SetTitle()   
    "如果文件类型为.sh文件   
    if  &filetype == 'cpp'  
       	call setline(1, "/*************************************************************************")   
        call append(line("."), "    > File Name: ".expand("%"))   
        call append(line(".")+1, "    > Author: bbezxcy")   
        call append(line(".")+2, "    > Mail:[e-mail] ")   
        call append(line(".")+3, "    > Created Time: ".strftime("%c"))   
        call append(line(".")+4, " ************************************************************************/")   
        call append(line(".")+5, "")  
    elseif  &filetype == 'c'  
       	call setline(1, "/*************************************************************************")   
        call append(line("."), "    > File Name: ".expand("%"))   
        call append(line(".")+1, "    > Author: bbezxcy")   
        call append(line(".")+2, "    > Mail:[e-mail] ")   
        call append(line(".")+3, "    > Created Time: ".strftime("%c"))   
        call append(line(".")+4, " ************************************************************************/")   
        call append(line(".")+5, "")  
    else 
    	call setline(1,"\#########################################################################")   
        call append(line("."), "\# File Name: ".expand("%"))   
        call append(line(".")+1, "\# Author: bbezxcy")   
        call append(line(".")+2, "\# mail: [e-mail]")   
        call append(line(".")+3, "\# Created Time: ".strftime("%c"))   
        call append(line(".")+4, "\#########################################################################")   
        if  &filetype == 'sh'  
            call append(line(".")+5, "\#!/bin/bash")   
        else
            call append(line(".")+5, "\#coding=utf-8")   
        endif
        call append(line(".")+6, "")   
    endif  
    if &filetype == 'cpp'  
        call append(line(".")+6, "#include <iostream>") 
        call append(line(".")+7, "#include <cstdio>")
        call append(line(".")+8, "using namespace std;")  
        call append(line(".")+9, "")  
    endif  
    if &filetype == 'c'  
    call append(line(".")+6, "#include <stdio.h>")  
    call append(line(".")+7, "")  
    endif  
    "新建文件后，自动定位到文件末尾  
    autocmd BufNewFile * normal G  
endfunc     
 
" 使用VIM的键盘 这样下面的配置才能有效
set nocompatible
set syntax=on   "设置语法高亮
set autoindent  "设置自动缩进
set cindent
set softtabstop=4   "统一缩进为 4
set shiftwidth=4
set tabstop=4   "设置tab 宽度为 4
set noexpandtab "不用空格代替制表符！ 便于删除
set number  "显示行号
 
"自动补全功能
:inoremap ( ()<ESC>i
:inoremap ) <c-r>=ClosePair(')')<CR>
:inoremap { {<CR>}<ESC>O
:inoremap } <c-r>=ClosePair('}')<CR>
:inoremap [ []<ESC>i
:inoremap ] <c-r>=ClosePair(']')<CR>
function! ClosePair(char)
    if getline('.')[col('.') - 1] == a:char
        return "\<Right>"
    else
        return a:char
    endif
endfunction
 
"C，C++ 按F5编译运行
map <F5> :call CompileRunGcc()<CR>
func! CompileRunGcc()
    exec "w"
    if &filetype == 'c'
        exec "!g++ % -o %<"
        exec "! ./%<"
    elseif &filetype == 'cpp'
        exec "!g++ % -o %<"
        exec "! ./%<"
    elseif &filetype == 'java'
        exec "!javac %"
        exec "!java %<"
    elseif &filetype == 'sh'
        :!./%
    elseif &filetype == 'py'
 	exec "!python %"
    endif
endfunc
 
"C,C++的 F8调试
map <F8> :call Rungdb()<CR>
func! Rungdb()
    exec "w"
    exec "!g++ % -g -o %<"
    exec "!gdb ./%<"
endfunc

 

评论

2 楼 tyfengyu 2015-03-19  

我刚刚更改的代码加上了标准差stdVal,故recoMat应该为：
recoMat = finalData * tfMat.T * stdVal+ meanVal

1 楼 tyfengyu 2015-03-19  

python的pca代码有2处错误：
1.finalData = rmmeanMat*tfMat 是点乘，numpy的矩阵乘应该用finalData =dot(rmmeanMat,tfMat)
2.tfMat =eigvec[0:lenth,:]  没有排序，该取列而不是行；


更改如下：
from numpy import * 
def pca(mat, length): 
    meanVal = mean(mat, axis = 0)
    stdVal = std(mat)
    rmmeanMat = (mat - meanVal)/stdVal
    covMat = cov(rmmeanMat,rowvar = 0) 
    eigval,eigvec =  linalg.eig(covMat)
    maxnum=argsort(-eigval, axis=0)#sort descend
    tfMat =eigvec[:,maxnum[0:length]]#top length
    finalData = dot(rmmeanMat,tfMat)#
    recoMat = finalData * tfMat.T + meanVal
    return finalData,recoMat