こんにちはゲストさん。会員登録(無料)して質問・回答してみよう!

締切り済みの質問

最小二乗法を応用した微分

C言語についてです。
規則性のある離散データを局所的(たとえば3点ごと)に最小二乗法を使って微分したいのですがどうすればいいでしょうか。
以下、データ全体に最小二乗法を適用するプログラムです。
#include<stdio.h>
#define n 60
main()
{
FILE *fin;
int i;
float x[n],y[n];
double a=0,b=0,sumx=0,sumy=0,sumxy=0,sumx2=0;
fin=fopen("data1.txt","r");

if(fin==NULL){
fprintf(stderr,"Not found\n");
}
for(i=0;!feof(fin);++i){
fscanf(fin,"%f %f",&x[i],&y[i]);
}

for(i=0;i<n;++i){
sumx+=x[i];
sumy+=y[i];
sumxy+=x[i]*y[i];
sumx2+=x[i]*x[i];
}
a=(n*sumxy-sumx*sumy)/(n*sumx2-sumx*sumx);
b=(sumx2*sumy-sumxy*sumx)/(n*sumx2-sumx*sumx);
printf("a=%f b=%f\n",a,b);
fclose(fin);
}

データは以下のようなものです。
0 0.5
0.1 0.8
0.2 1.4
・・・・・
回答お願いします。

投稿日時 - 2007-12-26 16:32:30

QNo.3629015

すぐに回答ほしいです

このQ&Aは役に立ちましたか?

0人が「このQ&Aが役に立った」と投票しています

回答(2)

ANo.2

方法(手順)ってアルゴリズムそのもののような気がしますが、、

例えば、
・ご質問中のプログラムを
 「データ入力」の部分と「最小二乗法の計算(と結果出力)」の部分に切り分ける。
・最小二乗法の計算部分を、関数にしてmainの外に出す。
 (このとき、引数として、データの開始点s、処理する個数 kを渡すとする。例えばsaishou_jijou(int s, int k) の形で)

という下準備をしておけば、
main(){
int i;
/*データ入力部分*/
.....
for(i=0;i<n-3;i++){
saishou_jijou(i,3);
}
みたいな形で処理できそうに思います。

(この方法だと、計算量が、O(元のデータ量*平均する区間のデータ量)程度になります。処理するデータ量が非常に多い場合には、最小二乗法の部分を工夫する必要がありますが、60点のデータから3点分の平均傾きを計算するくらいの処理量なら、特に工夫を加えなくても大丈夫かと思います。)
}
}

投稿日時 - 2007-12-28 06:38:56

補足

御回答ありがとうございます。
回答内容を参考に(ほぼそのままですが・・・)以下のようなプログラムを書いてみたのですが、うまくいきません。
ほぼ同じ値が出力されます。
どこかまちがっているのでしょうか。回答お願いします。
#include<stdio.h>
#define n 60
float x[n],y[n];
double a=0,b=0,sumx=0,sumy=0,sumxy=0,sumx2=0;
int calc(int p,int q)
{
int i;
for(i=0;i<n;++i){
sumx+=x[i];
sumy+=y[i];
sumxy+=x[i]*y[i];
sumx2+=x[i]*x[i];
}
a=(q*sumxy-sumx*sumy)/(q*sumx2-sumx*sumx);
b=(sumx2*sumy-sumxy*sumx)/(q*sumx2-sumx*sumx);
printf("a=%f b=%f\n",a,b);
}

main()
{
int i;
/*データ入力*/
・・・
for(i=0;!feof(fin);++i){
fscanf(fin,"%f %f",&x[i],&y[i]);
}
for(i=0;i<n-3;++i){
calc(i,3);
}
fclose(fin);
}

投稿日時 - 2008-01-04 17:02:46

ANo.1

単純には、
1.ある区間(3点?)で最小二乗法でa,bを求める。
2.aをその区間の(平均的な)微係数として採用する。
3.区間を1サンプル分移動して、1,2を行う
というの測定点の最初から最後まで行えば良さそうに思います。

「数値微分」のキーワードで検索すれば、他にもいろいろ参考となる資料が引っかかるかと思います。

投稿日時 - 2007-12-26 19:10:06

補足

説明不足でした。
方法はわかるのですが、そのアルゴリズムがわかりません。

投稿日時 - 2007-12-27 14:07:38

あなたにオススメの質問