什么网站都能打开的浏览器,网站详情页链接怎么做,asp大型网站开发,苏州网站开发公司招聘信息任务
iris数据集包含150条数据#xff0c;从iris.txt读取#xff0c;每条数据有4个属性值和一个标签#xff08;标签取值为0#xff0c;1#xff0c;2#xff09;。要求对这150个4维数据进行PCA#xff0c;可视化展示这些数据在前两个主方向上的分布#xff0c;其中不…任务
iris数据集包含150条数据从iris.txt读取每条数据有4个属性值和一个标签标签取值为012。要求对这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×nm1XXT的特征值选取最大两个特征值对应的特征向量构成矩阵 P 2 × n \mathbf{P}_{2\times n} P2×n则 Y 2 × m P X \mathbf{Y}_{2\times m}\mathbf{PX} Y2×mPX即PCA后的结果也就是把四维数据压缩为二维每个数据对应二维平面上的一个点。 对PCA的详解可以参考这篇文章关于PCA与奇异值分解的联系可以参考这篇文章关于如何用C求矩阵特征值Jacobi方法和特征向量及对矩阵进行奇异值分解可以参考这篇文章。
代码
C实现PCA
#includebits/stdc.h
using namespace std;// 读取鸢尾花数据集到一个二维数组中
vectorvectorlong double readIrisData(const string filename);// 读取第五列的值到一个向量中
vectorlong double readfifthValue(const string filename);// 从矩阵 A 非对角元中选择最大的元素并返回其位置
pairint, int chooseMax(vectorvectorlong double A);// 计算矩阵 A 的转置
vectorvectorlong double calAT(vectorvectorlong double A);// 计算矩阵 A 和其转置的乘积
vectorvectorlong double calAAT(vectorvectorlong double A);// 计算矩阵Q^T * A * Q的每个元素使用给定的参数 p, q, t, c, d
long double calculateElement(const vectorvectorlong 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
vectorvectorlong double calQTAQ(vectorvectorlong double A);// 判断Jacobi迭代方法是否满足结束条件
int judgeEnd(vectorvectorlong double A);// 计算矩阵 A 的特征值
vectorlong double calEigenValue(vectorvectorlong double A);// 对矩阵 A 进行列主元化成上三角
vectorvectorlong double Column_Elimination(vectorvectorlong double A);// 求解系数矩阵为上三角矩阵A的线性方程组
vectorlong double SolveUpperTriangle(vectorvectorlong double A, vectorlong double b);// 解系数矩阵为上三角矩阵 A 的线性方程组且A全为0的行数为 cnt
vectorvectorlong double solve(vectorvectorlong double A, int cnt);// 计算矩阵 A 的特征向量使用给定的特征值
vectorvectorlong double calEigenVector(vectorvectorlong double A, vectorlong double eigenValue);// 计算 Sigma 矩阵使用给定的特征值 x 和矩阵的行数 n1 和列数 n2
vectorvectorlong double calSigma(vectorlong double x,int n1, int n2);// 计算向量 x 的欧几里得范数
long double EuclideanNorm(vectorlong double x);// 对矩阵 A 进行归一化
vectorvectorlong double Normalization(vectorvectorlong double A);// 计算矩阵 A 和 B 的乘积
vectorvectorlong double multiplyMatrices(const vectorvectorlong double A, const vectorvectorlong double B);int main()
{vectorvectorlong 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;}vectorvectorlong double XT calAT(X);vectorvectorlong double XXT multiplyMatrices(X, XT);vectorvectorlong double Var(n1, vectorlong double(n1));for(int i 0; i n1; i)for(int j 0; j n1; j)Var[i][j] XXT[i][j] / n2;vectorlong 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] ;vectorlong double x1;x1.reserve(x.size());unique_copy(x.begin(), x.end(), back_inserter(x1));vectorvectorlong double EigenVector Normalization(calEigenVector(Var, x1));vectorvectorlong 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;}vectorvectorlong 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;
}// 读取鸢尾花数据集到一个二维数组中
vectorvectorlong double readIrisData(const string filename) {ifstream file(filename);vectorvectorlong double X;string line;while (getline(file, line)) {stringstream ss(line);vectorlong 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;
}// 读取第五列的值到一个向量中
vectorlong double readfifthValue(const string filename) {ifstream file(filename);vectorlong 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 中非对角元中绝对值最大的元素并返回其位置
pairint, int chooseMax(vectorvectorlong double A)
{long double max 0;pairint, 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 的转置
vectorvectorlong double calAT(vectorvectorlong double A)
{int n1 A.size();int n2 A[0].size();vectorvectorlong double AT(n2, vectorlong double(n1));for(int i 0; i n1; i)for(int j 0; j n2; j)AT[j][i] A[i][j];return AT;
}// 计算两个矩阵的乘积
vectorvectorlong double multiplyMatrices(const vectorvectorlong double A, const vectorvectorlong double B) {int n1 A.size();int n2 B[0].size();int n3 A[0].size();vectorvectorlong double result(n1, vectorlong 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 vectorvectorlong 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
vectorvectorlong double calQTAQ(vectorvectorlong double A)
{int n A.size();pairint, 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;vectorvectorlong double QTAQ(n, vectorlong 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(vectorvectorlong 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 的特征值
vectorlong double calEigenValue(vectorvectorlong double A)
{int n A.size();vectorlong double eigenValue(n);vectorvectorlong 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 进行列主元化成上三角
vectorvectorlong double Column_Elimination(vectorvectorlong double A)
{int n A.size();vectorvectorlong double Temp(n, vectorlong 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的线性方程组
vectorlong double SolveUpperTriangle(vectorvectorlong double A, vectorlong double b)
{int n A.size();vectorlong double x(n);vectorvectorlong double Temp(n, vectorlong double(n1));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
vectorvectorlong double solve(vectorvectorlong double A, int cnt)
{int n A.size();vectorvectorlong double x(cnt, vectorlong double(n));vectorvectorlong double Temp(n-cnt, vectorlong double(n-cnt));vectorlong 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];}vectorlong double res SolveUpperTriangle(Temp, Tempb);for(int j 0; j n - cnt; j)x[i][j] res[j];}return x;
}// 使用给定的特征值计算矩阵 A 的特征向量
vectorvectorlong double calEigenVector(vectorvectorlong double A, vectorlong double eigenValue)
{int n A.size();int num 0;vectorvectorlong double x(n, vectorlong double(n));vectorvectorlong double tempMartix(n, vectorlong double(n));vectorvectorlong double eigenVector(n, vectorlong 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];vectorvectorlong 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;}}vectorvectorlong 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 矩阵
vectorvectorlong double calSigma(vectorlong double x, int n1, int n2)
{vectorvectorlong double Sigma(n1, vectorlong double(n2));for(int i 0; i min(n1, n2); i)Sigma[i][i] sqrt(x[i]);return Sigma;
}// 计算向量 x 的欧几里得范数
long double EuclideanNorm(vectorlong double x)
{long double norm 0;for(int i 0; i x.size(); i)norm x[i] * x[i];return sqrt(norm);
}// 对矩阵 A 进行归一化
vectorvectorlong double Normalization(vectorvectorlong 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,, dtypefloat, encodingNone)return data[:, :4].T, data[:, 4]X, labels readIrisData(iris.txt)Var np.cov(X)
x, EigenVector eigh(Var)
x sorted(x, reverseTrue)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, colorcolors[i], markershapes[i], labellabel)plt.legend()#添加图例
plt.show()可视化结果
