#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define Nmax 100
#define Mmax 20
void main_part(char fnameR[],char fnameW[]);
void MATGJ(int n,double a[Mmax+1][Mmax+2]);
/*--------------------------------------------------------------------------*/

int main(int argc,char *argv[])
{
    char fnameR[50]; /*入力ファイル名*/
    char fnameW[50]; /*出力ファイル名*/

    strcpy(fnameR,argv[1]);
    strcpy(fnameW,argv[2]);

    main_part(fnameR,fnameW);

    printf("リターンキーを押してください");
    getchar();
    return 0;
}
/*---------------------------------------------------------------------------*/

void main_part(char fnameR[],char fnameW[])
{
    /*重回帰分析（multiple regression analysis）*/
    FILE *fin,*fout;
    char s[256];

    char scom[256];
    int ndata;     /* データ組数−１ */
    int mcol;      /* 説明変数の数−１ */
    int mm;        /* 説明変数の数−１＋１ */
    double RR;     /* 重相関係数 */
    double YEmean; /* 回帰推定値の平均 */
    double YDmean; /* 目的変数入力値の平均 */
    double std;    /* 残差の標準偏差 */
    double XD[Nmax+1][Mmax+1]; /* 説明変数入力値用配列 */
    double XT[Mmax+1][Nmax+1]; /* 説明変数入力値用配列の転置行列 */
    double YD[Nmax+1];         /* 目的変数入力値 */
    double a[Mmax+1][Mmax+2];  /* 連立方程式用配列 */
    double BB[Mmax+1];         /* 偏回帰係数 */
    double YE[Nmax+1];         /* 回帰推定値 */
    double RS[Nmax+1];         /* 目的変数入力値と回帰推定値の差分 */

    int i,j,k;
    double y1,y2,y3;

    /* データ読み込み */
    fin=fopen(fnameR,"r");
    fgets(s,sizeof s,fin);strcpy(scom,strtok(s,","));
    fgets(s,sizeof s,fin);mcol=atoi(strtok(s,","))-1;ndata=atoi(strtok(NULL,","))-1;
    mm=mcol+1;
    i=0;
    while(fgets(s,sizeof s,fin)!=NULL){
        XD[i][0]=1.0;
        XD[i][1]=atof(strtok(s,","));
        for(j=1;j<=mcol;j++){
            XD[i][j+1]=atof(strtok(NULL,","));
        }
        YD[i]=atof(strtok(NULL,","));
        i=i+1;
    }
    fclose(fin);

    /* 正規方程式作成 */
    for(i=0;i<=ndata;i++){
        for(j=0;j<=mm;j++){
            XT[j][i]=XD[i][j];
    }}
    for(i=0;i<=mm;i++){
        for(j=0;j<=mm;j++){
            a[i][j]=0.0;
            for(k=0;k<=ndata;k++){
                a[i][j]=a[i][j]+XT[i][k]*XD[k][j];
            }
    }}
    for(i=0;i<=mm;i++){
        a[i][mm+1]=0.0;
        for(k=0;k<=ndata;k++){
            a[i][mm+1]=a[i][mm+1]+XT[i][k]*YD[k];
        }
    }

    /* 連立一次方程式の解 */
    MATGJ(mm,a);
    for(i=0;i<=mm;i++){
        BB[i]=a[i][mm+1];
    }

    /* 重相関係数算定 */
    YDmean=0.0;YEmean=0.0;
    for(i=0;i<=ndata;i++){
        YE[i]=0.0;
        for(j=0;j<=mm;j++){
            YE[i]=YE[i]+XD[i][j]*BB[j];
        }
        YDmean=YDmean+YD[i];
        YEmean=YEmean+YE[i];
    }
    YDmean=YDmean/(double)(ndata+1);
    YEmean=YEmean/(double)(ndata+1);
    y1=0.0;y2=0.0;y3=0.0;
    for(i=0;i<=ndata;i++){
        y1=y1+(YD[i]-YDmean)*(YE[i]-YEmean);
        y2=y2+(YD[i]-YDmean)*(YD[i]-YDmean);
        y3=y3+(YE[i]-YEmean)*(YE[i]-YEmean);
    }
    RR=y1/sqrt(y2*y3);
    /* 残差の標準偏差 */
    y1=0.0;
    for(i=0;i<=ndata;i++){
        RS[i]=YD[i]-YE[i];
        y1=y1+RS[i];
    }
    y1=y1/(double)(ndata+1); /* 残差の平均 */
    y2=0.0;
    for(i=0;i<=ndata;i++){
        y2=y2+(RS[i]-y1)*(RS[i]-y1);
    }
    std=sqrt(y2/(double)(ndata)); /* 残差の標準偏差 */

    /* 分析結果の画面（ListBox）への書き出し */
    printf(scom);
    printf("y=B0+B1*x1+B2*x2+・・・+Bm*xm\n");
    for(i=0;i<=mm;i++){
        printf("B[%d]=%e\n",i,BB[i]);
    }
    printf("重相関係数：R=%e\n",RR);
    printf("データ数：n=%d\n",ndata+1);
    printf("残差平均：%e\n",y1);
    printf("残差標準偏差：%e\n",std);

    /*結果書き出し*/
    fout=fopen(fnameW,"w");
    fprintf(fout,"%s",scom);
    fprintf(fout,"y=B0+B1*x1+B2*x2+・・・+Bm*xm\n");
    fprintf(fout,"偏回帰係数\n");
    for(i=0;i<=mm;i++){
        fprintf(fout,"B[%d]=,%e\n",i,BB[i]);
    }
    fprintf(fout,"重相関係数：R=,%e\n",RR);
    fprintf(fout,"データ数：n=,%d\n",ndata+1);
    fprintf(fout,"残差平均：,%e\n",y1);
    fprintf(fout,"残差標準偏差：,%e\n",std);
    fprintf(fout,"目的変数入力値,回帰推定値\n");
    for(i=0;i<=ndata;i++){
        fprintf(fout,"%e,%e\n",YD[i],YE[i]);
    }
    fclose(fout);

}
/*--------------------------------------------------------------------------*/

void MATGJ(int n,double a[Mmax+1][Mmax+2])
{
    /* Gauss-Jordan法による連立一次方程式の解法 */
    int i,j,k,s;
    double p,d,max,dumy;

    for(k=0;k<=n;k++){
        max=0.0;
        s=k;
        for(j=k;j<=n;j++){
            if(fabs(a[j][k])>max){
                max=fabs(a[j][k]);
                s=j;
            }
        }
        for(j=0;j<=n+1;j++){
            dumy=a[k][j];
            a[k][j]=a[s][j];
            a[s][j]=dumy;
        }
        p=a[k][k];
        for(j=k;j<=n+1;j++){
            a[k][j]=a[k][j]/p;
        }
        for(i=0;i<=n;i++){
            if(i!=k){
                d=a[i][k];
                for(j=k;j<=n+1;j++){
                    a[i][j]=a[i][j]-d*a[k][j];
                }
            }
        }
    }
}
/*-------------------------------------------------------------------*/