回帰分析

回帰分析とは、データの関係をなるべくずれないように式で表すための手法である。

実装例

(x,y)の関係をy=a*x+bという式で表すとき、全ての点について(a*x+b-y)の二乗を求め、その和が最小になるようなa,bを求める。
この実装では、x,yの組を1行に一個ずつコンマ区切りで入力する。

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

int main(void) {
	/*入力用*/
	char inputbuf[100000];
	char* comma;
	char* comma2;
	/*入力チェック用*/
	char* endptr;
	/*統計用*/
	int datanum=0;
	int ignorenum=0;
	/*計算用*/
	double x,y;
	double abkeisuu;
	double ateisuu;
	double b2keisuu;
	double bkeisuu;
	double teisuu;
	double b;
	double a;
	double keisuu[6]={0,0,0,0,0,0};
	/*
	  0:a^2
	  1:b^2 
	  2:ab
	  3:a
	  4:b
	  5:定数
	*/
	/*データ入力*/
	while(fgets(inputbuf,sizeof(inputbuf),stdin)) {
		/*入力データの解析*/
		/*改行を消す*/
		comma=strchr(inputbuf,'\n');
		if(comma)*comma=0;
		/*一つ目のコンマ*/
		comma=strchr(inputbuf,',');
		if(comma==NULL) {
			/*コンマが一つもない*/
			ignorenum++;
			continue;
		}
		*comma=0;
		/*二つ目のコンマ*/
		comma2=strchr(comma+1,',');
		if(comma2)*comma2=0;/*いらないデータは消す*/
		/*変換とデータチェック*/
		x=strtod(inputbuf,&endptr);
		if(*endptr!=0) {
			ignorenum++;
			continue;
		}
		y=strtod(comma+1,&endptr);
		if(*endptr!=0) {
			ignorenum++;
			continue;
		}
		/*係数を加える*/
		keisuu[5]+=y*y;
		keisuu[0]+=x*x;
		keisuu[1]+=1;
		keisuu[3]-=2.0*x*y;
		keisuu[2]+=2.0*x;
		keisuu[4]-=2.0*y;
		/*データ数を増やす*/
		datanum++;
	}
	/*計算*/
	/*
	  a^2+b^2+ab+a+b+1
	  a^2+(b+1)a+b^2+b+1→2a+(b+1)=0→2a=-b-1→a=-b/2-1/2
	*/
	if(keisuu[0]==0) {
		puts("エラー：第一次0除算");
		return 1;
	}
	abkeisuu=-keisuu[2]/2.0/keisuu[0];
	ateisuu=-keisuu[3]/2.0/keisuu[0];
	b2keisuu=keisuu[0]*abkeisuu*abkeisuu+keisuu[1]+keisuu[2]*abkeisuu;
	bkeisuu=keisuu[0]*2.0*abkeisuu*ateisuu+keisuu[2]*ateisuu+keisuu[3]*abkeisuu+keisuu[4];
	teisuu=keisuu[0]*ateisuu*ateisuu+keisuu[3]*ateisuu+keisuu[5];
	if(b2keisuu==0) {
		puts("エラー：第二次0除算");
		return 1;
	}
	b=-bkeisuu/2.0/b2keisuu;
	a=abkeisuu*b+ateisuu;
	/*結果の出力*/
	printf("総行数：%d\n",datanum+ignorenum);
	printf("無視した行数：%d\n",ignorenum);
	printf("データの行数：%d\n",datanum);
	if(b>=0) {
		printf("y=%1.15gx+%1.15g\n",a,b);
	} else {
		printf("y=%1.15gx%1.15g\n",a,b);
	}
	return 0;
}

入力例(無効なデータは無視される)

時刻[s],速さ[m/s]
0.05,0.70
0.15,1.60
0.25,2.57
0.35,3.53
0.45,4.51

出力例

総行数：6
無視した行数：1
データの行数：5
y=9.55000000000001x+0.194499999999998

関連項目

• 最小自乗法
• 重回帰分析