PCA算法的主要步骤是:
让客户满意是我们工作的目标,不断超越客户的期望值来自于我们对这个行业的热爱。我们立志把好的技术通过有效、简单的方式提供给客户,将通过不懈努力成为客户在信息化领域值得信任、有价值的长期合作伙伴,公司提供的服务项目有:域名申请、网页空间、营销软件、网站建设、忻城网站维护、网站推广。
(1) 对向量X进行去中心化
(2) 计算向量X的协方差矩阵,自由度可以选择0或1
(3)计算协方差矩阵的特征值和特征向量
(4)选取最大的k个特征值及其特征向量
(5)用X与特征向量相乘
python实现:
from sklearn.datasets import load_iris
import numpy as np
def pca(X, k):
X = X - X.mean(axis=0)
X_cov = np.cov(X.T, ddof = 0)
eigenvalues, eigenvectors = eig(X_cov)
klarge_index = eigenvalues.argsort()[-k:][::-1]
k_eigenvectors = eigenvectors[klarge_index]
return np.dor(X, k_eigenvectors.T)
iris = load_iris()
X = iris.data
k = 2
X_pca = pca(X, k)
程序说明:
y = pca(mixedsig),程序中mixedsig为 n*T 阶混合数据矩阵,n为信号个数,T为采样点数, y为 m*T 阶主分量矩阵。
程序设计步骤:
1、去均值
2、计算协方差矩阵及其特征值和特征向量
3、计算协方差矩阵的特征值大于阈值的个数
4、降序排列特征值
5、去掉较小的特征值
6、去掉较大的特征值(一般没有这一步)
7、合并选择的特征值
8、选择相应的特征值和特征向量
9、计算白化矩阵
10、提取主分量
程序代码
%程序说明:y = pca(mixedsig),程序中mixedsig为 n*T 阶混合数据矩阵,n为信号个数,T为采样点数
% y为 m*T 阶主分量矩阵。
function y = pca(mixedsig)
if nargin == 0
error('You must supply the mixed data as input argument.');
end
if length(size(mixedsig))2
error('Input data can not have more than two dimensions. ');
end
if any(any(isnan(mixedsig)))
error('Input data contains NaN''s.');
end
%——————————————去均值————————————
meanValue = mean(mixedsig')';
mixedsig = mixedsig - meanValue * ones(1,size(meanValue,2));
[Dim,NumofSampl] = size(mixedsig);
oldDimension = Dim;
fprintf('Number of signals: %d\n',Dim);
fprintf('Number of samples: %d\n',NumofSampl);
fprintf('Calculate PCA...');
firstEig = 1;
lastEig = Dim;
covarianceMatrix = cov(mixedsig',1); %计算协方差矩阵
[E,D] = eig(covarianceMatrix); %计算协方差矩阵的特征值和特征向量
%———计算协方差矩阵的特征值大于阈值的个数lastEig———
rankTolerance = 1e-5;
maxLastEig = sum(diag(D)) rankTolerance;
lastEig = maxLastEig;
%——————————降序排列特征值——————————
eigenvalues = flipud(sort(diag(D)));
%—————————去掉较小的特征值——————————
if lastEig oldDimension
lowerLimitValue = (eigenvalues(lastEig) + eigenvalues(lastEig + 1))/2;
else
lowerLimitValue = eigenvalues(oldDimension) - 1;
end
lowerColumns = diag(D) lowerLimitValue;
%—————去掉较大的特征值(一般没有这一步)——————
if firstEig 1
higherLimitValue = (eigenvalues(firstEig - 1) + eigenvalues(firstEig))/2;
else
higherLimitValue = eigenvalues(1) + 1;
end
higherColumns = diag(D) higherLimitValue;
%—————————合并选择的特征值——————————
selectedColumns =lowerColumns higherColumns;
%—————————输出处理的结果信息—————————
fprintf('Selected[ %d ] dimensions.\n',sum(selectedColumns));
fprintf('Smallest remaining (non-zero) eigenvalue[ %g ]\n',eigenvalues(lastEig));
fprintf('Largest remaining (non-zero) eigenvalue[ %g ]\n',eigenvalues(firstEig));
fprintf('Sum of removed eigenvalue[ %g ]\n',sum(diag(D) .* (~selectedColumns)));
%———————选择相应的特征值和特征向量———————
E = selcol(E,selectedColumns);
D = selcol(selcol(D,selectedColumns)',selectedColumns);
%——————————计算白化矩阵———————————
whiteningMatrix = inv(sqrt(D)) * E';
dewhiteningMatrix = E * sqrt(D);
%——————————提取主分量————————————
y = whiteningMatrix * mixedsig;
%——————————行选择子程序———————————
function newMatrix = selcol(oldMatrix,maskVector)
if size(maskVector,1)~ = size(oldMatrix,2)
error('The mask vector and matrix are of uncompatible size.');
end
numTaken = 0;
for i = 1:size(maskVector,1)
if maskVector(i,1) == 1
takingMask(1,numTaken + 1) == i;
numTaken = numTaken + 1;
end
end
newMatrix = oldMatrix(:,takingMask);
本文来自CSDN博客,转载请标明出处:
/* --- 演示 STC 1T 系列单片机 用PCA功能实现16位定时器 --------*/
#include "reg51.h"
#include "intrins.h"
#define FOSC 18432000L
#define T100Hz (FOSC / 12 / 100)
typedef unsigned char BYTE;
typedef unsigned int WORD;
/*Declare SFR associated with the PCA */
sfr CCON = 0xD8; //PCA control register
sbit CCF0 = CCON^0; //PCA module-0 interrupt flag
sbit CCF1 = CCON^1; //PCA module-1 interrupt flag
sbit CR = CCON^6; //PCA timer run control bit
sbit CF = CCON^7; //PCA timer overflow flag
sfr CMOD = 0xD9; //PCA mode register
sfr CL = 0xE9; //PCA base timer LOW
sfr CH = 0xF9; //PCA base timer HIGH
sfr CCAPM0 = 0xDA; //PCA module-0 mode register
sfr CCAP0L = 0xEA; //PCA module-0 capture register LOW
sfr CCAP0H = 0xFA; //PCA module-0 capture register HIGH
sfr CCAPM1 = 0xDB; //PCA module-1 mode register
sfr CCAP1L = 0xEB; //PCA module-1 capture register LOW
sfr CCAP1H = 0xFB; //PCA module-1 capture register HIGH
sfr PCAPWM0 = 0xf2;
sfr PCAPWM1 = 0xf3;
sbit PCA_LED = P1^0; //PCA test LED
BYTE cnt;
WORD value;
void PCA_isr() interrupt 7 using 1
{
CCF0 = 0; //Clear interrupt flag
CCAP0L = value;
CCAP0H = value 8; //Update compare value
value += T100Hz;
if (cnt-- == 0)
{
cnt = 100; //Count 100 times
PCA_LED = !PCA_LED; //Flash once per second
}
}
void main()
{
CCON = 0; //Initial PCA control register
//PCA timer stop running
//Clear CF flag
//Clear all module interrupt flag
CL = 0; //Reset PCA base timer
CH = 0;
CMOD = 0x00; //Set PCA timer clock source as Fosc/12
//Disable PCA timer overflow interrupt
value = T100Hz;
CCAP0L = value;
CCAP0H = value 8; //Initial PCA module-0
value += T100Hz;
CCAPM0 = 0x49; //PCA module-0 work in 16-bit timer mode
//and enable PCA interrupt
CR = 1; //PCA timer start run
EA = 1;
cnt = 0;
while (1);
}
PCA的处理步骤:
1,均值化
2,求协方差矩阵(我知道的有两种方法,这是第一种,按部就班的求,第二种是:(A*A‘/(N-1)))
3,求协方差的特征值和特征向量
4,将特征值按照从大到小的顺序排序,选择其中最大的k个,然后将其对应的k个特征向量分别作为列向量组成特征向量矩阵
5,将样本点投影到选取的特征向量上
matlab实现源代码
%PCA算法,matlab实现
function F=pcad(A,n)%A是M*N
%测试实例A=[2.5,0.5,2.2,1.9,3.1,2.3,2,1,1.5,1.1;2.4,0.7,2.9,2.2,3.0,2.7,1.6,1.1,1.6,0.9]
%结果F=[0.8280,-1.7776,0.9922,0.2742,1.6758,0.9129,-0.0991,-1.1446,-0.4380,-1.2238]
%PCA第一步:均值化
X=A-repmat(mean(A,2),1,size(A,2))%去均值
%PCA第二步:求特征协方差矩阵
B=COV(X')%求协方差
%PCA第三步:求特征协方差矩阵的特征值和特征向量
[v,d]=eig(B)%求特征值和特征向量
%PCA第四步:将特征值按照从大到小的顺序排序
d1=diag(d);%取出对角矩阵,也就是把特征值提出来组成一个新的M*1的d1矩阵
[d2 index]=sort(d1); %特征值以升序排序 d2是排序后的结果 index是数排序以前的排名位置
cols=size(v,2);% 特征向量矩阵的列数
for i=1:cols %对特征向量做相反位置的调整 是个降序排列。这个过程把特征值和特征向量同时做相应的降序排列
vsort(:,i) = v(:,index(cols-i+1) ); % vsort 是一个M*col(注:col一般等于M)阶矩阵,保存的是按降序排列的特征向量,每一列构成一个特征向量
%vsort保存的是协方差矩阵降序后的特征向量,为M*M阶
dsort(i) = d1(index(cols-i+1)); % dsort 保存的是按降序排列的特征值,是一维行向量,1*M
end %完成降序排列
M=vsort(:,1:n)%提取主成分量
%PCA第五步:将样本点投影到选取的特征向量上
F=(X'*M)'%最终的投影
package PCA;
import java.util.ArrayList;
import java.util.Collections;
import java.util.List;
import org.jblas.ComplexDoubleMatrix;
import org.jblas.DoubleMatrix;
import org.jblas.Eigen;
public class PCA {
/**
* Reduce matrix dimension 减少矩阵维度
* @param source 源矩阵
* @param dimension 目标维度
* @return Target matrix 返回目标矩阵
*/
public static void main(String[] args){
DoubleMatrix d = new DoubleMatrix(new double[][]{{-1,-1,0,2,0},
{-2,0,0,1,1}});
DoubleMatrix result = PCA.dimensionReduction(d, 2);
System.out.println(result);
}
public static DoubleMatrix dimensionReduction(DoubleMatrix source, int dimension) {
//C=X*X^t/m 矩阵*矩阵^异或/列数
DoubleMatrix covMatrix = source.mmul(source.transpose()).div(source.columns);
ComplexDoubleMatrix eigVal = Eigen.eigenvalues(covMatrix);
ComplexDoubleMatrix[] eigVectorsVal = Eigen.eigenvectors(covMatrix);
ComplexDoubleMatrix eigVectors = eigVectorsVal[0];
//通过特征值将符号向量从大到小排序
ListPCABean beans = new ArrayListPCA.PCABean();
for (int i = 0; i eigVectors.columns; i++) {
beans.add(new PCABean(eigVal.get(i).real(), eigVectors.getColumn(i)));
}
Collections.sort(beans);
DoubleMatrix newVec = new DoubleMatrix(dimension, beans.get(0).vector.rows);
for (int i = 0; i dimension; i++) {
ComplexDoubleMatrix dm = beans.get(i).vector;
DoubleMatrix real = dm.getReal();
newVec.putRow(i, real);
}
return newVec.mmul(source);
}
static class PCABean implements ComparablePCABean {
double eigenValue;
ComplexDoubleMatrix vector;
public PCABean(double eigenValue, ComplexDoubleMatrix vector) {
super();
this.eigenValue = eigenValue;
this.vector = vector;
}
@Override
public int compareTo(PCABean o) {
return Double.compare(o.eigenValue, eigenValue);
}
@Override
public String toString() {
return "PCABean [eigenValue=" + eigenValue + ", vector=" + vector + "]";
}
}
}
org.jblas的jar包去百度下一个,我不知道怎么上传文件
private void lbl_CutImage_Paint(object sender, PaintEventArgs e)
{
int imgWidth = this.lbl_CutImage.Width - 4;
int imgHeight = this.lbl_CutImage.Height - 4;
if (imgWidth 1) { imgWidth = 1; }
if (imgHeight 1) { imgHeight = 1; }
// 创建缓存图像,先将要绘制的内容全部绘制到缓存中,最后再一次性绘制到 Label 上,
// 这样可以提高性能,并且可以防止屏幕闪烁的问题
Bitmap bmp_lbl = new Bitmap(this.lbl_CutImage.Width, this.lbl_CutImage.Height);
Graphics g = Graphics.FromImage(bmp_lbl);
// 将要截取的部分绘制到缓存
Rectangle destRect = new Rectangle(2, 2, imgWidth, imgHeight);
Point srcPoint = this.lbl_CutImage.Location;
srcPoint.Offset(2, 2);
Rectangle srcRect = new Rectangle(srcPoint, new System.Drawing.Size(imgWidth, imgHeight));
g.DrawImage(this.screenImage, destRect, srcRect, GraphicsUnit.Pixel);
SolidBrush brush = new SolidBrush(Color.FromArgb(10, 124, 202));
Pen pen = new Pen(brush, 1.0F);