liboctave for windows MSVC or linux g++ 2008
liboctave以外には,IT++などのライブラリがあり,(これはliboctaveと異なり,そもそもAPIとして設計されたものなので)ドキュメントが豊富です.この情報もメモしておきます.
以下に2001,2002年あたりの古い情報がありましたが,あまりにも古くなってきたので使用例を除いて削除しました.
Octave C++で数値計算例
固有値分解
EIG
00001 #include <iostream> 00002 #include <octave/config.h> 00003 #include <octave/Matrix.h> 00004 00005 using namespace std; 00006 00007 int main() 00008 { 00009 Matrix m(2, 2); 00010 m(0,0) = 3; m(0,1) = 2; 00011 m(1,0) = 2; m(1,1) = 0; 00012 cout << "Original Matrix" << endl << m << endl; 00013 00014 EIG eig(m); 00015 cout << "Eigen Vectors" << endl << eig.eigenvectors() << endl; 00016 cout << "Eigen Values" << endl << eig.eigenvalues() << endl; 00017 cout << "Recomposed Matrix" << endl << eig.eigenvectors() *ComplexMatrix(ComplexDiagMatrix(eig.eigenvalues())) *eig.eigenvectors().inverse() << endl; 00018 00019 return 0; 00020 }
特異値分解
00001 #include <iostream> 00002 #include <octave/config.h> 00003 #include <octave/Matrix.h> 00004 00005 using namespace std; 00006 00007 int main() 00008 { 00009 Matrix m(3, 2); 00010 m(0,0) = 3; m(0,1) = 2; 00011 m(1,0) = 2; m(1,1) = 0; 00012 m(2,0) = 4; m(2,1) = 2; 00013 cout << "Original Matrix" << endl << m << endl; 00014 00015 SVD svd(m); 00016 cout << "Left Singular Matrix" << endl << svd.left_singular_matrix() << endl; 00017 cout << "Singular Values" << endl << svd.singular_values() << endl; 00018 cout << "Right Singular Matrix" << endl << svd.right_singular_matrix() << endl; 00019 cout << "Recomposed Matrix" << endl << svd.left_singular_matrix() *Matrix(DiagMatrix(svd.singular_values())) *svd.right_singular_matrix() << endl; 00020 00021 return 0; 00022 }
Octave C++で多変量解析
平均と共分散・相関係数
00001 #include <iostream> 00002 #include <octave/config.h> 00003 #include <octave/Matrix.h> 00004 00005 #include "octaveio.h" 00006 00007 using namespace std; 00008 00009 int main(int argc, char **argv) 00010 { 00011 Matrix X; 00012 "X.dat" >> X; // 各列がサンプルデータ 00013 int dim = X.rows(); // ベクトルの次元数 00014 int nsample = X.cols(); // サンプルベクトル数 00015 cout << "X\n" << X << endl; 00016 00017 // 平均ベクトルの計算 00018 ColumnVector mean(dim); 00019 mean = (X.sum(1) / nsample).column(0); // X.sum(1)/nsample だとMatrix型のまま 00020 cout << "mean\n" << mean << endl; 00021 00022 // 共分散行列の計算 00023 Matrix Cov(dim, dim); 00024 Matrix Xd(dim, nsample); // 平均ベクトルを各列から引いた行列 00025 for(int i=0; i < nsample; i++){ 00026 Xd.insert(X.column(i) - mean, 0, i); 00027 } 00028 Cov = Xd*Xd.transpose() / nsample; // 不偏量を考えるなら nsample - 1 で割る 00029 cout << "Cov\n" << Cov << endl; 00030 00031 // 相関係数行列の計算 00032 Matrix Rate(dim, dim); 00033 Matrix Xn(dim, nsample); // 各行を正規化した行列 00034 for(int i=0; i < dim; i++){ 00035 Xn.insert(X.row(i)/sqrt(X.row(i)*X.row(i).transpose()), i, 0); 00036 } 00037 Rate = Xn*Xn.transpose(); 00038 cout << "Rate\n" << Rate << endl; 00039 00040 return 0; 00041 }
主成分分析
00001 #include <iostream> 00002 #include <octave/config.h> 00003 #include <octave/Matrix.h> 00004 00005 #include "octaveio.h" 00006 00007 using namespace std; 00008 00009 int main(int argc, char **argv) 00010 { 00011 Matrix X; 00012 "X.dat" >> X; // 各列がサンプルデータ 00013 int dim = X.rows(); // ベクトルの次元数 00014 int nsample = X.cols(); // サンプルベクトル数 00015 cout << "X\n" << X << endl; 00016 00017 ColumnVector mean(dim); // 平均ベクトル 00018 Matrix Xd(dim, nsample); // 平均ベクトルを各列から引いた行列 00019 Matrix Cov(dim, dim); // 共分散行列 00020 mean = (X.sum(1) / nsample).column(0); 00021 for(int i=0; i < nsample; i++){ 00022 Xd.insert(X.column(i) - mean, 0, i); 00023 } 00024 Cov = Xd*Xd.transpose() / nsample; 00025 00026 // 主成分分析(固有値分解を用いて) 00027 EIG eig(Cov); 00028 cout << "eigen vectors\n" << eig.eigenvectors() << endl; // 主成分 00029 cout << "eigen values\n" << eig.eigenvalues() << endl; // 固有値(昇順) 00030 00031 // 主成分分析(特異値分解を用いて) 00032 SVD svd(Xd); 00033 cout << "eigen vectors\n" << svd.left_singular_matrix() << endl; 00034 ColumnVector eigenvalues(dim); 00035 for(int i=0; i < dim; i++){ // Xdを直接分解したので共分散行列の固有値のサンプル数倍 00036 eigenvalues(i) = svd.singular_values()(i,i)*svd.singular_values()(i,i) / nsample; 00037 } 00038 cout << "eigen values\n" << eigenvalues << endl; // 固有値(降順) 00039 00040 return 0; 00041 }
正準相関分析
00001 #include <iostream> 00002 #include <octave/config.h> 00003 #include <octave/Matrix.h> 00004 00005 #include "octaveio.h" 00006 00007 using namespace std; 00008 00009 int main() 00010 { 00011 // 2組の2次元変数群(x_i, y_i (i=1,...,nsample))を読み込む 00012 // ただし平均が0であることを仮定する 00013 Matrix X; 00014 Matrix Y; 00015 "X.dat" >> X; 00016 "Y.dat" >> Y; 00017 00018 cout << "X\n" << X << endl; 00019 cout << "Y\n" << Y << endl; 00020 00021 int dim = X.rows(); 00022 int nsample = X.cols(); 00023 00024 // それぞれの変数群を白色化するために,主成分分析を行う 00025 Matrix Ux, Uy; 00026 DiagMatrix Sx_inv, Sy_inv; 00027 EIG eig_X(X*X.transpose()/nsample); 00028 EIG eig_Y(Y*Y.transpose()/nsample); 00029 Ux = real(eig_X.eigenvectors()); 00030 Uy = real(eig_Y.eigenvectors()); 00031 ColumnVector Sx = real(eig_X.eigenvalues()); 00032 ColumnVector Sy = real(eig_Y.eigenvalues()); 00033 Sx_inv.resize(dim,dim); 00034 Sy_inv.resize(dim,dim); 00035 for(int i=0; i < dim; i++){ 00036 Sx_inv(i,i) = 1.0/sqrt(Sx(i)); 00037 Sy_inv(i,i) = 1.0/sqrt(Sy(i)); 00038 } 00039 00040 // 白色化する行列 (XX'/nsample)^(-1/2), (YY'/nsample)^(-1/2) 00041 Matrix Wx = Ux*(Matrix)Sx_inv*Ux.transpose(); 00042 Matrix Wy = Uy*(Matrix)Sy_inv*Uy.transpose(); 00043 00044 // 白色化された変量群 00045 Matrix Xw = Wx * X; 00046 Matrix Yw = Wy * Y; 00047 00048 // 白色化された変量群間の相互相関を特異値分解する 00049 SVD svd_R(Yw*Xw.transpose()/nsample, SVD::economy); 00050 // 係数行列 00051 Matrix A = svd_R.right_singular_matrix().transpose() * Wx; 00052 Matrix B = svd_R.left_singular_matrix().transpose() * Wy; 00053 cerr << "A\n" << A << endl; 00054 // 正準相関係数 00055 DiagMatrix Sr = svd_R.singular_values(); 00056 cerr << "Sr\n" << Sr << endl; 00057 00058 //-------------------------------------------------- 00059 // 検証 00060 //-------------------------------------------------- 00061 // 正準変量群 00062 Matrix Xc = A * X; 00063 Matrix Yc = B * Y; 00064 cout << "Xc\n" << Xc << endl; 00065 cout << "Yc\n" << Yc << endl; 00066 // 正準変量群間の相互相関は,正準相関係数が対角に並んだ行列 00067 cout << "Xc*Yc'/nsample\n" << Xc*Yc.transpose()/nsample << endl; 00068 00069 return 0; 00070 }