#include <stdio.h>
#include <math.h>
#define Nmax 10
void KAIKI(int ndata, double datax[Nmax+1], double datay[Nmax+1],
           double *aa, double *bb, double *rr);

/*---------------------------------------------------------------------------*/

int main()
{
    int ndata = 7;
    double datax[Nmax];
    double datay[Nmax];
    double aa,bb,rr;

    datax[0]=1.0;datay[0]=3.0;
    datax[1]=2.0;datay[1]=5.0;
    datax[2]=3.0;datay[2]=7.0;
    datax[3]=4.0;datay[3]=9.0;
    datax[4]=5.0;datay[4]=11.0;
    datax[5]=6.0;datay[5]=13.0;
    datax[6]=7.0;datay[6]=15.0;
    datax[7]=8.0;datay[7]=17.0;

    KAIKI(ndata, datax, datay, &aa, &bb, &rr);

    printf("y=a*x+b\n");
    printf("a=%e\n",aa);
    printf("b=%e\n",bb);
    printf("r=%e\n",rr);

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

void KAIKI(int ndata, double datax[Nmax], double datay[Nmax],
           double *aa, double *bb, double *rr)
{
    int i;
    double x1,y1,x2,xy,xm,ym,c1,c2,c3;
            
    /* 回帰式：y=aa*x+bb */
    x1 = 0.0 ; y1 = 0.0 ; x2 = 0.0 ; xy = 0.0;
    for(i=0;i<=ndata;i++){
        x1 = x1 + datax[i];
        y1 = y1 + datay[i];
        x2 = x2 + datax[i] * datax[i];
        xy = xy + datax[i] * datay[i];
    }
    xm = x1 / (double)(ndata + 1) ; ym = y1 / (double)(ndata + 1);
    *aa = ((double)(ndata + 1) * xy - x1 * y1) / ((double)(ndata + 1) * x2 - x1 * x1);
    *bb = (x2 * y1 - x1 * xy) / ((double)(ndata + 1) * x2 - x1 * x1);
    c1 = 0.0 ; c2 = 0.0 ; c3 = 0.0;
    for(i=0;i<=ndata;i++){
        c1 = c1 + (datax[i] - xm) * (datay[i] - ym);
        c2 = c2 + (datax[i] - xm) * (datax[i] - xm);
        c3 = c3 + (datay[i] - ym) * (datay[i] - ym);
    }
    *rr = c1 / sqrt(c2 * c3);
}
/*---------------------------------------------------------------------------*/