プログラムなど個人的メモ

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

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

エクセルを用いた最小二乗法 : http://homepage1.nifty.com/gfk/square_ren.htm
逆行列の求め方 : http://thira.plavox.info/blog/2008/06/_c.html

プログラム

//行列の作成
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;
	    }
	}    	
    }
}

とりあえずテストしてみて、エクセルのLINEST関数の結果と照らし合わせるなどして比較。いけそうだと思ったがうまくいかない場合がでてきた。その場合とは小数点の値を使う場合。たとえば、x = 0.1、0.2、0.3 を使用しての４次元の多項式y=a * x^4 + b * x^3 + c * x^2 + d * x^1 + e の各a,b,c,d,eを求めるために、各項を2乗、３乗としていったが、デバッグで値を見ると、後ろに妙な端数がある。
これは、浮動小数点型を使用して計算を行うとどうしても誤差がでるらしい。
http://dobon.net/vb/dotnet/beginner/floatingpointerror.html

多少誤差が出ようが、似たような計算結果になってほしいところだが、実際に上記のプログラムを使用して算出された値を使い、使用したx値を与えてやっても同じようなy値が算出できなかった。

なので、結局エクセルのLinEstクラスを使用するようにした。
使用するには、HDF5DotNetの時と同じく、参照の追加を行わなければならない。追加するものは、COMにある"Microsoft Office xx(xxはofficeのバージョン)Object Library"。これを追加すると、ソリューションエクスプローラーの中の参照設定の中に"Microsoft.Office.Core"と"Microsoft.Office.Interop.Excel"が追加される。

これでエクセルの編集などができるようになったが、あくまで目的はワークシート関数の"LinEst"を使用すること。

とりあえず、簡単に以下のラッパークラスを作成した。

static public  class CoefficientsClass
{
    static Excel.Application oXL = new Excel.Application( );
    static Excel.WorksheetFunction funk = oXL.WorksheetFunction;
    
    public static Array GetCoefficients(Array y, Array x, bool a, bool b)
    {			
	System.Object coefficients = funk.LinEst( y, x, a, b );
	Array coeffArray = (Array)coefficients;
	return coeffArray;
    }
}

説明するとstaticでエクセル起動とワークシート関数をクラスに置いておく。
データを持ってきてLinEst関数を使用する。ただ単にラップしただけ。
あとはこれを使うだけ

double[ , ] x = new double[ dim , xyValues.Count ];
double[] y = new double[ xyValues.Count ];
for( int i = 0 ; i < xyValues.Count ; i++ ) {
    for( int j = 0 ; j < dim ; j++ ) {
	x[ j, i ] = Math.Round(Math.Pow(xyValues[i].Key,j + 1),
		15,MidpointRounding.AwayFromZero);
    }
    y[ i ] = xyValues[ i ].Value;
}
 
Array coeffArray = CoefficientsClass.GetCoefficients( y, x ,true,false);
List<double> ret = new List<double>();
for(int i = 1 ; i < coeffArray.Length + 1 ; i++){
    ret.Add( (double)coeffArray.GetValue( i ) );
}

xを配列単体で行うとy=ax+bの係数と切片が得られる。xを上記のように二次配列にして、ほしい次元まで乗数倍したx値を用意すれば、多項近似式も得ることができる。