HDF5DotNetとは別でプログラムの話。

プログラムで、多項式近似の係数を出力する機能がほしくなったので、適当にネットで転がっている情報をもとに作成。参考ページとできたのがこれ

エクセルを用いた最小二乗法 : [URL]
逆行列の求め方 : [URL]

プログラム

//行列の作成
dimWithIntercept = 5;
double[ , ] matrix = new double[ dimWithIntercept, dimWithIntercept];
double[ , ] copy = new double[ dimWithIntercept, dimWithIntercept ];
double[ , ] revMatrix = new double[ dimWithIntercept, dimWithIntercept ];
double[] vector = new double[dimWithIntercept];
List<double> ret = new List<double>( );
int digit = 15;
MidpointRounding rounding = MidpointRounding.AwayFromZero;
for( int i = 0 ; i < xyValues.Count ; i++ ) {
    //       行列の作成
    for( int j = 0 ; j < dimWithIntercept ; j++ ) {
        for( int k = 0 ; k < dimWithIntercept ; k++ ) {
            matrix[ j, k ] = matrix[ j, k ] + Math.Pow( xyValues[ i ].Key, j );
            copy[ j, k ] = matrix[ j, k ];
        }
//ベクトルの作成
    vector[ j ] = vector[ j ] + Math.Pow(xyValues[ i ].Key, j ) ;
    }                           
}
//   逆行列の作成
//   まずは単位行列を作る
for( int i = 0 ; i < dimWithIntercept ; i++ ) {
    for( int j = 0 ; j < dimWithIntercept ; j++ ) {
        revMatrix[ i,j ] = (i == j) ? 1.0 : 0.0;
    }
}
//      吐き出し法
double buf;
for( int i = 0 ; i < dimWithIntercept ; i++ ) {
    //       単位化するために、対になる値で割る
    buf = copy[ i, i ];
    //       iの行を単位化するためにbufで割る
for( int j = 0 ; j < dimWithIntercept ; j++ ) {
        copy[ i, j ] = copy[ i, j ] / buf;
        revMatrix[ i, j ] = revMatrix[ i, j ] / buf;
    }
    for( int j = 0 ; j < dimWithIntercept ; j++ ) {
    //        iと同じ行では計算しない
        if( i != j ) {
            buf = copy[ j, i ];
            //n列目の値を0にする(n,nの場所はスルー)
            for( int k = 0 ; k < dimWithIntercept ; k++ ) {
                double a = copy[ i, k ] * buf;
                copy[ j, k ] = copy[ j, k ] - a;
            a = revMatrix[ i, k ] * buf;
                revMatrix[ j, k ] =revMatrix[ j, k ] - a;