以鸢尾花数据为例,用C语言实现PCA算法,并求出4个主成分的贡献率。
时间: 2023-12-03 13:43:05 浏览: 42
PCA算法是一种常用的降维算法,可以将高维数据降到低维,保留数据的主要信息。以下是用C语言实现PCA算法的代码:
```c
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define N 150 // 样本数
#define M 4 // 特征数
// 求平均值
void mean(double data[N][M], double mean[M])
{
int i, j;
for (j = 0; j < M; j++) {
mean[j] = 0;
for (i = 0; i < N; i++) {
mean[j] += data[i][j];
}
mean[j] /= N;
}
}
// 求协方差矩阵
void cov(double data[N][M], double cov[M][M], double mean[M])
{
int i, j, k;
for (i = 0; i < M; i++) {
for (j = 0; j < M; j++) {
cov[i][j] = 0;
for (k = 0; k < N; k++) {
cov[i][j] += (data[k][i] - mean[i]) * (data[k][j] - mean[j]);
}
cov[i][j] /= (N - 1);
}
}
}
// 求特征值和特征向量
void eig(double cov[M][M], double eigval[M], double eigvec[M][M])
{
int i, j;
double eps = 1e-8;
double tmp[M][M];
for (i = 0; i < M; i++) {
for (j = 0; j < M; j++) {
tmp[i][j] = cov[i][j];
}
}
// 对称矩阵迭代法
while (1) {
double max = 0;
int p, q;
for (i = 0; i < M; i++) {
for (j = i + 1; j < M; j++) {
double aii = tmp[i][i];
double ajj = tmp[j][j];
double aij = tmp[i][j];
double bij = 2 * aij / (aii - ajj);
double t = bij > 0 ? -bij + sqrt(1 + bij * bij) : -bij - sqrt(1 + bij * bij);
double c = 1 / sqrt(1 + t * t);
double s = t * c;
double aik, ajk;
for (k = 0; k < M; k++) {
aik = tmp[i][k];
ajk = tmp[j][k];
tmp[i][k] = c * aik - s * ajk;
tmp[j][k] = s * aik + c * ajk;
}
for (k = 0; k < M; k++) {
aik = tmp[k][i];
ajk = tmp[k][j];
tmp[k][i] = c * aik - s * ajk;
tmp[k][j] = s * aik + c * ajk;
}
}
}
// 判断是否收敛
for (i = 0; i < M; i++) {
for (j = i + 1; j < M; j++) {
if (fabs(tmp[i][j]) > max) {
max = fabs(tmp[i][j]);
p = i;
q = j;
}
}
}
if (max < eps) {
break;
}
}
// 特征向量矩阵
for (i = 0; i < M; i++) {
for (j = 0; j < M; j++) {
eigvec[j][i] = tmp[j][i];
}
}
// 特征值向量
for (i = 0; i < M; i++) {
eigval[i] = tmp[i][i];
}
}
int main()
{
int i, j;
double data[N][M];
// 读取数据
FILE *fp = fopen("iris.data", "r");
if (!fp) {
printf("Open file failed!\n");
return -1;
}
for (i = 0; i < N; i++) {
fscanf(fp, "%lf,%lf,%lf,%lf", &data[i][0], &data[i][1], &data[i][2], &data[i][3]);
}
fclose(fp);
// 求平均值
double mean[M];
mean(data, mean);
// 求协方差矩阵
double cov[M][M];
cov(data, cov, mean);
// 求特征值和特征向量
double eigval[M];
double eigvec[M][M];
eig(cov, eigval, eigvec);
// 输出特征值和特征向量
printf("Eigenvalues:\n");
for (i = 0; i < M; i++) {
printf("%lf ", eigval[i]);
}
printf("\n");
printf("Eigenvectors:\n");
for (i = 0; i < M; i++) {
for (j = 0; j < M; j++) {
printf("%lf ", eigvec[i][j]);
}
printf("\n");
}
// 求贡献率
double sum = 0;
for (i = 0; i < M; i++) {
sum += eigval[i];
}
double ratio[M];
for (i = 0; i < M; i++) {
ratio[i] = eigval[i] / sum * 100;
}
// 输出贡献率
printf("Ratios:\n");
for (i = 0; i < M; i++) {
printf("%lf%% ", ratio[i]);
}
printf("\n");
return 0;
}
```
这个代码读取了一个名为`iris.data`的文件,文件中存储了150个鸢尾花的4个特征(花萼长度、花萼宽度、花瓣长度、花瓣宽度),按行存储。代码中先求出每个特征的平均值,然后计算协方差矩阵,最后通过对称矩阵迭代法求出特征值和特征向量。输出特征值和特征向量后,再根据特征值求出每个主成分的贡献率。