计算方法实验5：对鸢尾花数据集进行主成分分析(PCA)并可视化

2024-04-15 18:28
文章标签 分析 进行 数据 可视化 实验 成分 计算方法 鸢尾花 pca

本文主要是介绍计算方法实验5：对鸢尾花数据集进行主成分分析(PCA)并可视化，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

任务

iris数据集包含150条数据，从iris.txt读取，每条数据有4个属性值和一个标签（标签取值为0，1，2）。要求对这150个4维数据进行PCA，可视化展示这些数据在前两个主方向上的分布，其中不同标签的数据需用不同的颜色或形状加以区分。

算法

m个n维数据向量去中心化后（各向量的每个维度减去这个维度在所有向量上均值），按列排列构成矩阵 X n × m \mathbf{X}_{n\times m} Xn×m​，计算协方差矩阵 V a r n × n = 1 m X X T \mathbf{Var}_{n\times n}= \frac{1}{m}\mathbf{XXT} Varn×n​=m1​XXT的特征值，选取最大两个特征值对应的特征向量构成矩阵 P 2 × n \mathbf{P}_{2\times n} P2×n​，则 Y 2 × m = P X \mathbf{Y}_{2\times m}=\mathbf{PX} Y2×m​=PX即PCA后的结果，也就是把四维数据压缩为二维，每个数据对应二维平面上的一个点。
对PCA的详解，可以参考这篇文章；关于PCA与奇异值分解的联系，可以参考这篇文章；关于如何用C++求矩阵特征值（Jacobi方法）和特征向量及对矩阵进行奇异值分解，可以参考这篇文章。

代码

C++实现PCA

#include<bits/stdc++.h>
using namespace std;// 读取鸢尾花数据集到一个二维数组中
vector<vector<long double>> readIrisData(const string& filename);// 读取第五列的值到一个向量中
vector<long double> readfifthValue(const string& filename);// 从矩阵 A 非对角元中选择最大的元素，并返回其位置
pair<int, int> chooseMax(vector<vector<long double>> A);// 计算矩阵 A 的转置
vector<vector<long double>> calAT(vector<vector<long double>> A);// 计算矩阵 A 和其转置的乘积
vector<vector<long double>> calAAT(vector<vector<long double>> A);// 计算矩阵Q^T * A * Q的每个元素，使用给定的参数 p, q, t, c, d
long double calculateElement(const vector<vector<long double>> A, int i, int j, long double p, long double q, long double t, long double c, long double d);// 计算矩阵 Q^T * A * Q
vector<vector<long double>> calQTAQ(vector<vector<long double>> A);// 判断Jacobi迭代方法是否满足结束条件
int judgeEnd(vector<vector<long double>> A);// 计算矩阵 A 的特征值
vector<long double> calEigenValue(vector<vector<long double>> A);// 对矩阵 A 进行列主元化成上三角
vector<vector<long double>> Column_Elimination(vector<vector<long double>> A);// 求解系数矩阵为上三角矩阵A的线性方程组
vector<long double> SolveUpperTriangle(vector<vector<long double>> A, vector<long double> b);// 解系数矩阵为上三角矩阵 A 的线性方程组，且A全为0的行数为 cnt
vector<vector<long double>> solve(vector<vector<long double>> A, int cnt);// 计算矩阵 A 的特征向量，使用给定的特征值
vector<vector<long double>> calEigenVector(vector<vector<long double>> A, vector<long double> eigenValue);// 计算 Sigma 矩阵，使用给定的特征值 x 和矩阵的行数 n1 和列数 n2
vector<vector<long double>> calSigma(vector<long double> x,int n1, int n2);// 计算向量 x 的欧几里得范数
long double EuclideanNorm(vector<long double> x);// 对矩阵 A 进行归一化
vector<vector<long double>> Normalization(vector<vector<long double>> A);// 计算矩阵 A 和 B 的乘积
vector<vector<long double>> multiplyMatrices(const vector<vector<long double>> A, const vector<vector<long double>> B);int main()
{vector<vector<long double>> X = calAT(readIrisData("iris.txt"));int n1 = X.size();int n2 = X[0].size();long double sum = 0;for(int i = 0; i < n1; i++){long double sum = 0;for(int j = 0; j < n2; j++)sum += X[i][j];long double avg = sum / n2;for(int j = 0; j < n2; j++)X[i][j] -= avg;}cout << "X: " << endl;for(int i = 0; i < n1; i++){for(int j = 0; j < n2; j++)cout << X[i][j] << " ";cout << endl;}vector<vector<long double>> XT = calAT(X);vector<vector<long double>> XXT = multiplyMatrices(X, XT);vector<vector<long double>> Var(n1, vector<long double>(n1));for(int i = 0; i < n1; i++)for(int j = 0; j < n1; j++)Var[i][j] = XXT[i][j] / n2;vector<long double> x =calEigenValue(Var);sort(x.begin(), x.end());reverse(x.begin(), x.end());cout<<"特征值："<<endl;for(int i = 0; i < n1; i++)cout << x[i] << " ";vector<long double> x1;x1.reserve(x.size());unique_copy(x.begin(), x.end(), back_inserter(x1));vector<vector<long double>> EigenVector = Normalization(calEigenVector(Var, x1));vector<vector<long double>> P(EigenVector.begin(), next(EigenVector.begin(), 2));cout << "P: " << endl;for(int i = 0; i < 2; i++){for(int j = 0; j < n1; j++)cout << P[i][j] << " ";cout << endl;}vector<vector<long double>> Y = multiplyMatrices(P, X);cout << "Y: " << endl;for(int i = 0; i < 2; i++){for(int j = 0; j < n2; j++)cout << Y[i][j] << " ";cout << endl;}return 0;
}// 读取鸢尾花数据集到一个二维数组中
vector<vector<long double>> readIrisData(const string& filename) {ifstream file(filename);vector<vector<long double>> X;string line;while (getline(file, line)) {stringstream ss(line);vector<long double> row;string value;int counter = 0;while (getline(ss, value, ',') && counter < 4) {row.push_back(stod(value));counter++;}X.push_back(row);}return X;
}// 读取第五列的值到一个向量中
vector<long double> readfifthValue(const string& filename) {ifstream file(filename);vector<long double> fifthValues;string line;while (getline(file, line)) {stringstream ss(line);string value;int counter = 0;while (getline(ss, value, ',') && counter < 4) {counter++;}if (counter == 4) { long double fifthValue = stold(value);fifthValues.push_back(fifthValue);}}return fifthValues;
}// 找到矩阵 A 中非对角元中绝对值最大的元素，并返回其位置
pair<int, int> chooseMax(vector<vector<long double>> A)
{long double max = 0;pair<int, int> maxPos;int n = A.size();for(int i = 0; i < n; i++)for(int j = 0; j < n; j++)if(i != j && fabsl(A[i][j]) > max){max = fabsl(A[i][j]);maxPos = make_pair(i, j);}return maxPos;
}// 计算矩阵 A 的转置
vector<vector<long double>> calAT(vector<vector<long double>> A)
{int n1 = A.size();int n2 = A[0].size();vector<vector<long double>> AT(n2, vector<long double>(n1));for(int i = 0; i < n1; i++)for(int j = 0; j < n2; j++)AT[j][i] = A[i][j];return AT;
}// 计算两个矩阵的乘积
vector<vector<long double>> multiplyMatrices(const vector<vector<long double>> A, const vector<vector<long double>> B) {int n1 = A.size();int n2 = B[0].size();int n3 = A[0].size();vector<vector<long double>> result(n1, vector<long double>(n2, 0.0));for(int i = 0; i < n1; i++) {for(int j = 0; j < n2; j++) {for(int k = 0; k < n3; k++) {result[i][j] += A[i][k] * B[k][j];}}}return result;
}// 计算矩阵Q^T * A * Q的每个元素，使用给定的参数 p, q, t, c, d
long double calculateElement(const vector<vector<long double>> A, int i, int j, long double p, long double q, long double t, long double c, long double d) {if (i == p && j == p)return A[p][p] - t * A[p][q];else if (i == q && j == q)return A[q][q] + t * A[p][q];else if ((i == p && j == q) || (i == q && j == p))return 0;else if (i != q && i != p && (j == p || j == q))return (j == p ? c : d) * A[p][i] - (j == p ? d : (-c)) * A[q][i];else if ((i == p || i == q) && j != q && j != p)return (i == p ? c : d) * A[p][j] - (i == p ? d : (-c)) * A[q][j];elsereturn A[i][j];
}// 计算矩阵 Q^T * A * Q
vector<vector<long double>> calQTAQ(vector<vector<long double>> A)
{int n = A.size();pair<int, int> maxPos = chooseMax(A);int row = maxPos.first;int col = maxPos.second;long double s = (A[col][col] - A[row][row]) / (2 * A[row][col]);long double t = 0;if(s == 0)t = 1;else if(abs(-s + sqrt(1 + s * s)) <= abs(-s - sqrt(1 + s * s)))t = -s + sqrt(1 + s * s);elset = -s - sqrt(1 + s * s);long double c = 1 / sqrt(1 + t * t);long double d = t * c;vector<vector<long double>> QTAQ(n, vector<long double>(n));for(int i = 0; i < n; i++)for(int j = 0; j < n; j++)QTAQ[i][j] = calculateElement(A, i, j, row, col, t, c, d);return QTAQ;
}// 判断Jacobi迭代方法是否满足结束条件
int judgeEnd(vector<vector<long double>> A)
{int i, j;int n = A.size();for(i = 0; i < n; i++)for(j = 0; j < n; j++)if(i != j && fabsl(A[i][j]) >= 1e-6)return 0;if(i == n && j == n) return 1;
}// 计算矩阵 A 的特征值
vector<long double> calEigenValue(vector<vector<long double>> A)
{int n = A.size();vector<long double> eigenValue(n);vector<vector<long double>> QTAQ= calQTAQ(A);int i, j;while(!judgeEnd(QTAQ))QTAQ = calQTAQ(QTAQ);for(i = 0; i < n; i++)eigenValue[i] =QTAQ[i][i];return eigenValue;
}// 对矩阵 A 进行列主元化成上三角
vector<vector<long double>> Column_Elimination(vector<vector<long double>> A)
{int n = A.size();vector<vector<long double>> Temp(n, vector<long double>(n));for(int i = 0; i < n; i++)for(int j = 0; j < n; j++)Temp[i][j] = A[i][j];for(int col = 0; col < n; col++){long double maxnum = abs(Temp[col][col]);int maxrow = col;for(int row = col + 1; row < n; row++){if(abs(Temp[row][col]) > maxnum){maxnum = abs(Temp[row][col]);maxrow = row;}}swap(Temp[col], Temp[maxrow]);for(int row = col + 1; row < n; row++){long double res = Temp[row][col] / Temp[col][col];for(int loc = col; loc < n; loc++)Temp[row][loc] -= Temp[col][loc] * res; }}return Temp;
}// 求解系数矩阵为上三角矩阵A的线性方程组
vector<long double> SolveUpperTriangle(vector<vector<long double>> A, vector<long double> b)
{int n = A.size();vector<long double> x(n);vector<vector<long double>> Temp(n, vector<long double>(n+1));for(int i = 0; i < n; i++)for(int j = 0; j < n; j++)Temp[i][j] = A[i][j];for(int i = 0; i < n; i++)Temp[i][n] = b[i];for(int row = n-1; row >= 0; row--){for(int col = row + 1; col < n; col++){Temp[row][n] -= Temp[col][n] * Temp[row][col] / Temp[col][col];Temp[row][col] = 0;}Temp[row][n] /= Temp[row][row];Temp[row][row] = 1;}for(int i = 0; i < n; i++)x[i] = Temp[i][n];return x;
}// 解系数矩阵为上三角矩阵 A 的线性方程组，且A全为0的行数为 cnt
vector<vector<long double>> solve(vector<vector<long double>> A, int cnt)
{int n = A.size();vector<vector<long double>> x(cnt, vector<long double>(n));vector<vector<long double>> Temp(n-cnt, vector<long double>(n-cnt));vector<long double> Tempb(n-cnt);for(int i = 0; i < cnt; i++){for(int j = n - 1; j >= n - cnt; j--){if(j >= n - i)x[i][j] = 0;elsex[i][j] = 1;}}for(int i = 0; i < n - cnt; i++)for(int j = 0; j < n - cnt; j++)Temp[i][j] = A[i][j];for(int i = 0; i < cnt; i++){for(int j = n - cnt - 1; j >=  0; j--){Tempb[j] = 0;for(int k = 0; k < cnt; k++)Tempb[j] -= A[j][n- cnt + k] * x[i][n- cnt + k];}vector<long double> res = SolveUpperTriangle(Temp, Tempb);for(int j = 0; j < n - cnt; j++)x[i][j] = res[j];}return x;
}// 使用给定的特征值计算矩阵 A 的特征向量
vector<vector<long double>> calEigenVector(vector<vector<long double>> A, vector<long double> eigenValue)
{int n = A.size();int num = 0;vector<vector<long double>> x(n, vector<long double>(n));vector<vector<long double>> tempMartix(n, vector<long double>(n));vector<vector<long double>> eigenVector(n, vector<long double>(n));for(int k = 0; k < n; k++){for(int i = 0; i < n; i++)for(int j = 0; j < n; j++)i == j ? tempMartix[i][j] = A[i][j] - eigenValue[k] : tempMartix[i][j] = A[i][j];vector<vector<long double>> B = Column_Elimination(tempMartix);int cnt = 0;//记录消元后全为0的行数for(int i = 0; i < n; i++){for(int j = 0; j < n; j++){if(fabsl(B[i][j]) > 1e-7)break;else if(j == n - 1)cnt++;}}vector<vector<long double>> result = solve(B, cnt);for(int i = 0; i < cnt; i++)copy(result[i].begin(), result[i].end(), x[num + i].begin());num += cnt;}return x;
}// 使用给定的特征值 x 和矩阵的行数 n1 和列数 n2，计算 Sigma 矩阵
vector<vector<long double>> calSigma(vector<long double> x, int n1, int n2)
{vector<vector<long double>> Sigma(n1, vector<long double>(n2));for(int i = 0; i < min(n1, n2); i++)Sigma[i][i] = sqrt(x[i]);return Sigma;
}// 计算向量 x 的欧几里得范数
long double EuclideanNorm(vector<long double> x)
{long double norm = 0;for(int i = 0; i < x.size(); i++)norm += x[i] * x[i];return sqrt(norm);
}// 对矩阵 A 进行归一化
vector<vector<long double>> Normalization(vector<vector<long double>> A)
{int rows = A.size();for(int i = 0; i < rows; i++){long double norm = EuclideanNorm(A[i]);int cols = A[i].size();for(int j = 0; j < cols; j++)A[i][j] /= norm;}return A;
}

python实现PCA并将结果可视化

import numpy as np
from scipy.linalg import eigh
import matplotlib.pyplot as pltdef readIrisData(filename):data = np.genfromtxt(filename, delimiter=',', dtype='float', encoding=None)return data[:, :4].T, data[:, 4]X, labels = readIrisData("iris.txt")Var = np.cov(X)
x, EigenVector = eigh(Var)
x = sorted(x, reverse=True)P = EigenVector[:, -2:].T
Y = np.dot(P, X)plt.figure()
label_set = set(labels)
colors = ['r', 'g', 'b']
shapes = ['o', 's', '^']for i, label in enumerate(label_set):#enumerate函数返回每个标签及其索引x = [Y[0, j] for j in range(Y.shape[1]) if labels[j] == label]y = [Y[1, j] for j in range(Y.shape[1]) if labels[j] == label]plt.scatter(x, y, color=colors[i], marker=shapes[i], label=label)plt.legend()#添加图例
plt.show()

可视化结果

在这里插入图片描述

这篇关于计算方法实验5：对鸢尾花数据集进行主成分分析(PCA)并可视化的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！


http://www.chinasem.cn/article/906576。

相关文章

如何使用celery进行异步处理和定时任务(django)

如何使用celery进行异步处理和定时任务(django)

《如何使用celery进行异步处理和定时任务(django)》文章介绍了Celery的基本概念、安装方法、如何使用Celery进行异步任务处理以及如何设置定时任务,通过Celery,可以在Web应用中... 目录一、celery的作用二、安装celery三、使用celery 异步执行任务四、使用celery
阅读更多...
详谈redis跟数据库的数据同步问题

详谈redis跟数据库的数据同步问题

《详谈redis跟数据库的数据同步问题》文章讨论了在Redis和数据库数据一致性问题上的解决方案,主要比较了先更新Redis缓存再更新数据库和先更新数据库再更新Redis缓存两种方案,文章指出,删除R... 目录一、Redis 数据库数据一致性的解决方案1.1、更新Redis缓存、删除Redis缓存的区别二
阅读更多...
Redis事务与数据持久化方式

Redis事务与数据持久化方式

《Redis事务与数据持久化方式》该文档主要介绍了Redis事务和持久化机制,事务通过将多个命令打包执行,而持久化则通过快照（RDB）和追加式文件（AOF）两种方式将内存数据保存到磁盘,以防止数据丢失... 目录一、Redis 事务1.1 事务本质1.2 数据库事务与redis事务1.2.1 数据库事务1.
阅读更多...
Oracle Expdp按条件导出指定表数据的方法实例

Oracle Expdp按条件导出指定表数据的方法实例

《OracleExpdp按条件导出指定表数据的方法实例》：本文主要介绍Oracle的expdp数据泵方式导出特定机构和时间范围的数据,并通过parfile文件进行条件限制和配置,文中通过代码介绍... 目录1.场景描述 2.方案分析3.实验验证 3.1 parfile文件3.2 expdp命令导出4.总结
阅读更多...
Redis主从复制实现原理分析

Redis主从复制实现原理分析

《Redis主从复制实现原理分析》Redis主从复制通过Sync和CommandPropagate阶段实现数据同步,2.8版本后引入Psync指令,根据复制偏移量进行全量或部分同步,优化了数据传输效率... 目录Redis主DodMIK从复制实现原理实现原理Psync: 2.8版本后总结Redis主从复制实
阅读更多...
更改docker默认数据目录的方法步骤

更改docker默认数据目录的方法步骤

《更改docker默认数据目录的方法步骤》本文主要介绍了更改docker默认数据目录的方法步骤,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友们下面随着小编来一... 目录1.查看docker是否存在并停止该服务2.挂载镜像并安装rsync便于备份3.取消挂载备份和迁
阅读更多...
不删数据还能合并磁盘? 让电脑C盘D盘合并并保留数据的技巧

不删数据还能合并磁盘? 让电脑C盘D盘合并并保留数据的技巧

《不删数据还能合并磁盘?让电脑C盘D盘合并并保留数据的技巧》在Windows操作系统中，合并C盘和D盘是一个相对复杂的任务，尤其是当你不希望删除其中的数据时，幸运的是，有几种方法可以实现这一目标且在... 在电脑生产时，制造商常为C盘分配较小的磁盘空间，以确保软件在运行过程中不会出现磁盘空间不足的问题。但在
阅读更多...
锐捷和腾达哪个好? 两个品牌路由器对比分析

锐捷和腾达哪个好? 两个品牌路由器对比分析

《锐捷和腾达哪个好?两个品牌路由器对比分析》在选择路由器时，Tenda和锐捷都是备受关注的品牌，各自有独特的产品特点和市场定位，选择哪个品牌的路由器更合适，实际上取决于你的具体需求和使用场景，我们从... 在选购路由器时，锐捷和腾达都是市场上备受关注的品牌，但它们的定位和特点却有所不同。锐捷更偏向企业级和专
阅读更多...
SpringBoot使用minio进行文件管理的流程步骤

SpringBoot使用minio进行文件管理的流程步骤

《SpringBoot使用minio进行文件管理的流程步骤》MinIO是一个高性能的对象存储系统,兼容AmazonS3API,该软件设计用于处理非结构化数据,如图片、视频、日志文件以及备份数据等,本文... 目录一、拉取minio镜像二、创建配置文件和上传文件的目录三、启动容器四、浏览器登录 minio五、
阅读更多...
Java如何接收并解析HL7协议数据

Java如何接收并解析HL7协议数据

《Java如何接收并解析HL7协议数据》文章主要介绍了HL7协议及其在医疗行业中的应用,详细描述了如何配置环境、接收和解析数据,以及与前端进行交互的实现方法,文章还分享了使用7Edit工具进行调试的经... 目录一、前言二、正文1、环境配置2、数据接收：HL7Monitor3、数据解析：HL7Busines
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2