回帰分析
2013年9月5日 (木) 06:39時点におけるimported>MikeCATによる版 (作成)
回帰分析とは、データの関係をなるべくずれないように式で表すための手法である。
実装例
(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