FC2ブログ

最小二乗法のN次多項式近似

昨日の最小二乗法の1次多項式近似を最小二乗法のN次多項式近似に一般化したものです。


#define MAXDIMN 6
static void gauss(int dim, double a[MAXDIMN][MAXDIMN+1],double *outx)
{
int i,j,k,l;

for (i = 0; i < dim; i++) {
double m = 0;
int pivot = i;
for (l = i; l < dim; l++) {
// max
if (fabs(a[l][i]) > m) {
m = fabs(a[l][i]);
pivot = l;
}
}

if (pivot != i) {
// 行の入れ替え
for (j = 0; j < dim+1; j++) {
double t = a[i][j];
a[i][j] = a[pivot][j];
a[pivot][j] = t;
}
}
}

for (k = 0; k < dim; k++) {
double p = a[k][k];// 対角要素
a[k][k] = 1; // 対角要素は1になる
for (j = k+1; j < dim+1; j++) {
a[k][j] /= p;
}
for (i = k+1; i < dim; i++) {
double q = a[i][k];
for (j = k+1; j < dim+1; j++) {
a[i][j] -= q*a[k][j];
}
a[i][k] = 0; // 0となる
}
}

//
double x[MAXDIMN];
for (i = dim-1; i >= 0; i--) {
x[i] = a[i][dim];
for (j = dim-1; j>i; j--) {
x[i] -= a[i][j]*x[j];
}
}

for(i = 0; i < dim; i++) {
outx[i] = x[i];
}
}

void LeastSquaresMethod(double *x,double *y, int n, int dim, double *coefficient)
{
int i,j,k;
double A[MAXDIMN][MAXDIMN+1];

for (i = 0; i < dim; i++) {
for (j = 0; j < dim+1; j++) {
A[i][j] = 0.0;
}
}
for (i = 0; i < dim; i++) {
for (j = 0; j < dim; j++) {
for ( k = 0; k < n; k++) {
A[i][j] += pow(x[k],i+j);
}
}
}
for (i = 0; i < dim; i++) {
for (k = 0; k < n; k++) {
A[i][dim] +=pow(x[k],i)*y[k];
}
}

gauss(dim, A, coefficient);
}


3次での呼び出しは
dim=4で
double coefficient[4]を渡せばOK
返却データは、順に0次から3次の係数

double *x = new double[n];
double *y = new double[n];
double coefficient[4];
for (int i = 0; i < n; i++) {
x[i] = xxx;
y[i] = yyy;
}
LeastSquaresMethod(x, y, n, 4, coefficient);

スポンサーサイト

コメントの投稿

非公開コメント

プロフィール

chattera.net

Author:chattera.net
フリーでソフト開発の仕事をしています。
自作したこのチャットシステムは、何か世の中に役に立つこともあるかと思いますので、お呼びがかかるとうれしいですね。

連絡先は chatteranet@gmail.com です。

また、ソフト開発関連のお誘いもうけたまわっております。

最新記事
最新コメント
最新トラックバック
月別アーカイブ
カテゴリ
検索フォーム
RSSリンクの表示
リンク
ブロとも申請フォーム

この人とブロともになる

QRコード
QR