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

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

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;
    }
}

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 ) );
}