//-----------------------------------------------------------------
// corrmat.c:
//    相関行列の計算
//-----------------------------------------------------------------
#include <stdio.h>
#include <math.h>

//-----------------------------------------------------------------
double* calc_corrmat(int n,int samp_n,double* data);
double  calc_corr(int n,double* x,double* y);
double  calc_mean(int n,double* y);
void    disp_matrix(int m,int n,double* R);
//-----------------------------------------------------------------
static char fmt[] = "%.3lf "; // 表示フォーマット

int main()
{
    // <サンプルデータ>
    //工場の品質管理部門で, ある化学薬品の精製工程における呈色度 x1は,
    //その薬品に含まれるある微量不純物質の量 x2 と含湿度 x3 によって
    //近似的に線型の影響を受けると考えられていた. このことを現在の改良
    //された反応条件下で検討するため, 15 ロットからサンプルを無作為に
    //抽出し, (x1, x2, x3) を測定した. データはこの観測値を定数倍した
    //ものである. 
    double sample_data[45] = {
        53,33,27,63,23,43,63,43,53,23,13,40,30,43,40,         // x1
        50,60,10,53,13,57,7,20,10,37,7,6,40,53,23,            // x2
        130,77,90,110,73,97,137,113,130,100,87,127,133,93,123 // x3
    };
    double* R;
    // disp original data
    printf("**sample data:\n");
    disp_matrix(3,15,sample_data);
    // calc corr matrix
    R = calc_corrmat(3,15,sample_data);
    
    // disp corr matrix
    printf("**corr matrix:\n");
    disp_matrix(3,3,R);

    free(R);
    return 0;
}
// calc corr matrix
double* calc_corrmat(int n,int samp_n,double* data)
{
    int i,j;
    double *R = (double*)calloc(n*n,sizeof(double));
    // calc corr factor
    for(i=0;i<n;i++) for(j=i+1;j<n;j++){
        R[j+i*n] = calc_corr(samp_n,data+i*samp_n,data+j*samp_n);
        R[i+j*n] = R[j+i*n];
    }
    // diagonal factor
    for(i=0;i<n;i++) R[i+i*n] = 1.0;
    return R;
}

// calc corr
double calc_corr(int n,double* x,double* y)
{
    int i;
    double mx=calc_mean(n,x);
    double my=calc_mean(n,y);
    double xx=0,yy=0,xy=0,corr=0;
    for(i=0;i<n;i++){
        xx+=(x[i]-mx)*(x[i]-mx);
        yy+=(y[i]-my)*(y[i]-my);
        xy+=(x[i]-mx)*(y[i]-my);
    }
    xx = sqrt(xx); yy = sqrt(yy); corr = xy / (xx*yy);
    return corr;
}
// calc mean
double calc_mean(int n,double* a)
{
    int i;
    double mean=0;
    for(i=0;i<n;i++) mean += a[i];
    mean/=n;
    return mean;
}
// disp matrix
void disp_matrix(int m,int n,double* R)
{
    int i,j;
    for(i=0;i<m;i++){
        for(j=0;j<n;j++) printf(fmt,R[j+i*n]);
        printf("\n");
    }
}