Wednesday, July 12, 2023

PCAのコード

 毎回書き直すので、基本的なPCAを計算するコードを載せておこう。前も載せた気がするけど、ま、いいや。

c++でライブラリにEigenを使っている。

実験として、2次元10個のデータのmatのPCAを計算する。

とりあえず、ランダムのデータにする。

Eigen::MatrixXd mat(10, 2);

for (int i = 0; i < mat.rows(); i++) {

for (int k = 0; k < mat.cols(); k++) {

mat(i, k) = (double)(rand() % 100) / 10. + 50.*(double)k;

}

}


本編は以下で、PCAを使って2次元を1次元にプロジェクションする。

Eigen::MatrixXd aligned = mat.rowwise() - mat.colwise().mean();

Eigen::JacobiSVD<Eigen::MatrixXd> svd(aligned, Eigen::ComputeThinV);

Eigen::MatrixXd W = svd.matrixV().leftCols(1);

//これが主成分方向の内積の集合のベクトル

Eigen::MatrixXd dot_on_plane = aligned * W;

Eigen::MatrixXd projected_one = dot_on_plane * W.transpose();

Eigen::MatrixXd projected = projected_one.rowwise() + mat.colwise().mean();


No comments: