/*　このファイルは、「Cによる統計データの解析入門」で説明したプログラムのファイルである。
program1.txtには1-7章のプログラムが、program2に8-13章のプログラムがテキスト形式で収録されている。
これらのファイルは、通常の運用に関してはなんら問題のないことを確認しているが、運用はすべて自己責任で行うものとする。
運用の結果、万一損害が発生したとしても著者および出版社はいかなる責任も負わない。
これらのファイルに収録されたプログラムおよびソースコードの著作権は著者に属する。
ただし、読者が個人的に使用する場合および非営利団体において教育目的に使われる場合に関しては複写、改変、その一部の流用は自由である。
商用利用に関しては著者の同意が必要である。特に、著者の文書による同意なしに収録されたプログラムおよびソースコードを使って有料の著作物や製品の作成、
サービスの供与などを行うことは固く禁じる。詳細は出版社に連絡頂きたい。*/


/* 8章のプログラム　*/

/* Program 8-1(データ分析のプログラム) */

#include <stdio.h>
#include <math.h>
#include <string.h>
#define N 10

void inputdata(int num[ ], double x[ ], double y[ ], double z[ ], int n);
void calc_index(double x[ ], int n, char c[ ]);
double sum(double x[ ], int n);
double mean(double x[ ], int n); 
double variance(double x[ ], int n);
double stndv(double x[ ], int n);
double skewness(double x[ ],int n);
double kurtosis(double x[ ],int n);
void sort1(double x[ ], double y[ ], int n);
void quartile(double x[ ], double y[ ],int n, double *min, double *max,
			  double *q1, double *med, double *q3); 
double covar(double x[ ], double y[ ], int n);
double correl(double x[ ], double y[ ], int n);

void main( )
{
 int num[N];
 double x[N],y[N],z[N]; 
 char c[80];
 double covxy,rxy;
 
inputdata(num,x,y,z,N);
   strcpy(c,"Hight");
   calc_index(x,N,c);
   strcpy(c,"Weight");
   calc_index(y,N,c);
 covxy=covar(x,y,10);
 rxy=correl(x,y,10); 
 strcpy(c,"Hight and Weight");
 printf("\nCovariance and Correlation Coefficient: Variables: %s\n",c);
 printf("Covariance:%10.5f Correlaton Coefficient: %10.5f\n",covxy,rxy);
}

void calc_index(double x[ ], int n, char c[ ])
{double xsum, xmean,xvar,xstndv,xskew,xkur,xmin,xmax,xq1,
        xmed,xq3,xsort[N];
   xsum=sum(x,n);
   xmean=mean(x,n);
   xvar=variance(x,n);
   xstndv=stndv(x,n);
   xskew=skewness(x,n);
   xkur=kurtosis(x,n);
   quartile(x,xsort,n,&xmin,&xmax,&xq1,&xmed,&xq3);
printf("\nStatiatics of Variable: %s \n", c);  
printf("Sum:%10.3f Mean:%10.3f Variance:%10.3f Std Dev:%10.3f Skewness:%10.5f, Kurtosis:%10.5f\n",
	 xsum,xmean,xvar,xstndv,xskew,xkur);
printf("Minimum:%10.3f Maximum:%10.3f 1st Qrt%10.3f Mediam%10.3f, 3rd Qrt%10.3f\n\n",
	 xmin,xmax,xq1,xmed,xq3);
 }

void inputdata(int num[ ], double x[ ], double y[ ], double z[ ], int n)
{
int a1,i;
double a2,a3,a4;
FILE *fdata;
 fdata=fopen("c:\\cprog\\data1.txt","r");
   for(i=0; i<=9; i++)
	{
	 fscanf(fdata,"%d%lf%lf%lf",&a1,&a2,&a3,&a4);
	num[i]=a1;
	x[i]=a2;
	y[i]=a3;
	z[i]=a4;
	}
  fclose(fdata);
 }

double sum(double x[ ], int n)
{int i;
 double a;
   a=0;
   for(i=0; i<=n-1; i++)a=a+x[i];
return a;
}

double mean(double x[ ], int n)
{double a,b;
    a=sum(x,n);
    b=a/n;
 return b;
}

double variance(double x[ ], int n)
{double a,b,c;
int i;
   a=mean(x,n);
   b=0;
   for(i=0; i<=n-1; i++)b=b+pow(x[i]-a,2);
   c=b/(n-1);
return c;
}

double stndv(double x[ ], int n)
{double a,b;
   a=variance(x,n);
   b=sqrt(a);
return b;
}

double skewness(double x[ ], int n)
{double a,b,c,d;
int i;
   a=mean(x,n);
   b=stndv(x,n);
   c=0;
   for(i=0; i<=n-1; i++)c=c+pow(x[i]-a,3);
   d=c*n/(pow(b,3)*(n-1)*(n-2));
return d;
}

double kurtosis(double x[ ], int n)
{double a,b,c,d;
int i;
   a=mean(x,n);
   b=stndv(x,n);
   c=0;
   for(i=0; i<=n-1; i++)c=c+pow(x[i]-a,4);
   d=c*n*(n+1)/(pow(b,4)*(n-1)*(n-2)*(n-3));
   d=d-3.0*pow((n-1),2)/((n-2)*(n-3));
return d;
}

void sort1(double x[ ], double y[ ], int n)
{int change, i;
 double a;
    for(i=0; i<=n-1; i++) y[i]=x[i];
 do{
    change = 0;
 for(i=0; i<= n-2; i++)
    if(y[i]>y[i+1])
	{change=1;
	 a=y[i];
	 y[i]=y[i+1];
	 y[i+1]=a;}}
 while(change == 1);
}
void quartile(double x[ ], double y[ ], int n, double *min, double *max,
			  double *q1, double *med, double *q3) 
{int i,j1,j2;
   sort1(x,y,n);
   *min=y[0];
   *max=y[n-1];
     i=n%2;
     if(i!=0) *med=y[(n+1)/2-1];
     else *med=(y[n/2-1]+y[n/2])/2;
  i=n%4;
  j1=n/4;
  j2=3*n/4;
    if(i!=0){*q1=y[j1];
		*q3=y[j2];}
    else {*q1=(y[j1-1]+y[j1])/2;
       *q3=(y[j2-1]+y[j2])/2;}
}

double covar(double x[ ], double y[ ], int n)
{double xa,ya,a,b;
 int i;
   xa=mean(x,n);
   ya=mean(y,n);
   a=0;
   for(i=0; i<=9; i++)a=a+(x[i]-xa)*(y[i]-ya);
   b=a/(n-1);
 return b;
}

double correl(double x[ ], double y[ ], int n)
{double covxy,xstndv,ystndv,a;
   covxy=covar(x,y,n);
   xstndv=stndv(x,n);
   ystndv=stndv(y,n);
   a=covxy/(xstndv*ystndv);
return a;
}


/* 9章のプログラム　*/

/* Program 9-1 (二項乱数)　*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
int birnd(int n, float p);
int bern(float p);

void main( )
{int i,k,n,numrnd;
 float p;
 FILE *fout;
 fout=fopen("binrnd.txt","w");
   numrnd=200;
   n=5;
   p=0.5;
 for(i=1; i<=numrnd; i++)
 {k=birnd(n,p);
 fprintf(fout,"%5d\n",k);}
    fclose(fout);
}

int birnd(int n, float p)
{int i,k;
  k=0;
  for(i=1; i<=n; i++)k=k+bern(p);
  return k;
}

int bern(float p)
{int j;
double a;
a=rand( )/32768.0;
 if(a<p)j=1;
  else j=0;
  return j;
}


/* Program 9-2 (正規乱数) */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
double normrnd(void);
double norminv(double a);

void main( )
{int i,numrnd;
 double a;
 FILE *fout;
 fout=fopen("normrnd.txt","w");
   numrnd=1000;
 for(i=1; i<=numrnd; i++)
 {a=normrnd( );
 fprintf(fout,"%10f\n",a);}
    fclose(fout);
}

double normrnd(void)
{double a,b;
  a=rand()/32768.0;
  if(a==0) a=1/32768.0;
  b=norminv(a);
  return b;
}

double norminv(double u)
{double a,b,w,z;
  z=-log(4*u*(1-u));
  a=z*(2.0611786-5.7262204/(z+11.640595));
  w=sqrt(a);
  if(u<0.5)b=-w;
  else b=w;
return b;
}


/* Program 9-3 (指数乱数)　*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
double exprnd(double a);

void main( )
{int i,numrnd;
 double a;
 FILE *fout;
 fout=fopen("exprnd.txt","w");
   numrnd=200;
 for(i=1; i<=numrnd; i++)
 {a=exprnd(1);
 fprintf(fout,"%10f\n",a);}
    fclose(fout);
}

double exprnd(double a)
{double b,c;
  b=rand( )/32768.0;
  if(b==0) b=1/32768.0;
  c=-(1/a)*log(b);
  return c;
}


/* Program 9-4 (ポアソン乱数) */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
int pornd(double a);
double exprnd(double a);

void main( )
{int a,i,numrnd;
 FILE *fout;
 fout=fopen("pornd.txt","w");
 srand((unsigned)time(NULL));
   numrnd=200;
 for(i=1; i<=numrnd; i++)
 {a=pornd(3);
 fprintf(fout,"%5d\n",a);}
    fclose(fout);
}

int pornd(double t)
{double t1;
 int k;
 k=0;
  t1=exprnd(1);
  while(t1<t){
  k=k+1;
  t1=t1+exprnd(1);
  }
    return k;
}

double exprnd(double a)
{double b,c;
  b=rand( )/32768.0;
  if(b==0) b=1/32768.0;
  c=-(1/a)*log(b);
  return c;
}


/* Program 9-5 (大数の法則)　*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
double mean1(int n);

void main( )
{int i,n,numrnd;
 double a;
 FILE *fout;
 fout=fopen("LLN.txt","w");
   numrnd=200;
   n=5;
 for(i=1; i<=numrnd; i++)
 {a=mean1(n);
 fprintf(fout,"%10f\n",a);
 printf("%10f\n",a);}
    fclose(fout);
}

double mean1(int n)
{ double a,b,c;
  int i;
  a=0;
    for(i=1; i<=n; i++)
	{b=rand( )/32768.0;
    a=a+b;}
  c=a/n;
  return c;
}


/* Program 9-6 (中心極限定理) */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
double mean1(int n);
double mean2(int n);

void main( )
{int i,n,numrnd;
 double a;
 FILE *fout;
 fout=fopen("clt.txt","w");
   numrnd=200;
   n=5;
 for(i=1; i<=numrnd; i++)
 {a=mean2(n);
 fprintf(fout,"%10f\n",a);
 printf("%10f\n",a);}
    fclose(fout);
}

double mean1(int n)
{ double a,b,c;
  int i;
  a=0;
    for(i=1; i<=n; i++)
	{b=rand()/32768.0;
    a=a+b;}
  c=a/n;
  return c;
}

double mean2(int n)
{ double a,b;
  a=mean1(n);
  b=sqrt(12.0*n)*(a-0.5);  
  return b;
}


/* 10章のプログラム　*/

/* Program 10-1(行列を計算するプログラム) */

#include <stdio.h>
#include <math.h>
#define M 3
#define N 3
void inputdata(FILE *fdata, double x[][N],int m, int n);
void matsum(double x[ ][N], double y[ ][N], double z[ ][N], int m, int n);
void matprod(double x[ ][M], double y[ ][N], double z[ ][N], int m, int k, int n);
void mattrans(double x[ ][M], double y[ ][N], int m, int n);
void out_res(double a[ ][N], int m, int n);

void main( )
{ double a[3][N], b[3][N], c[2][N],sum1[3][N],prod1[3][N],trans1[3][N];
FILE *fdata;
  fdata=fopen("mat1.txt","r");
inputdata(fdata,a,3,2);
inputdata(fdata,b,3,2);
inputdata(fdata,c,2,3);
  matsum(a,b,sum1,3,2);
  matprod(a,c,prod1,3,2,3);
  mattrans(a,trans1,3,2);
	printf("  A+B \n");
    out_res(sum1,3,2);
	printf("  A*B \n");
    out_res(prod1,3,3);
	printf(" A' \n");
    out_res(trans1,2,3);
   fclose(fdata);
}

void inputdata(FILE *fdata,double x[ ][N],int m, int n)
{int i,j;
double d;
  for(i=0; i<= m-1; i++)
	  for(j=0; j<= n-1; j++)
	  {fscanf(fdata, "%lf", &d);
	  x[i][j]=d;
	  }
}

void matsum(double x[ ][N], double y[ ][N], double z[ ][N], int m, int n)
{int i,j;
 for(i=0; i<=m-1; i++)
	 for(j=0; j<=n-1; j++)
		 z[i][j]=x[i][j]+y[i][j];
}

void matprod(double x[][M], double y[ ][N], double z[ ][N], int m, int k, int n)
{int i,j,r;
 double a;
for(i=0; i<= m-1; i++)
for(j=0; j<= n-1; j++)
{   a=0;
	for(r=0; r<=k-1; r++) a=a+x[i][r]*y[r][j];
    z[i][j]=a;
}
}

void mattrans(double x[ ][M], double y[ ][N], int m, int n)
{int i,j;
for(i=0; i<=m-1; i++)
for(j=0; j<=n-1; j++)
y[j][i]=x[i][j];
}

void out_res(double x[ ][N], int m, int n)
{int i,j;
   for(i=0; i<= m-1; i++)
		{for(j=0; j<= n-1; j++)printf(" %10.3f ", x[i][j]);
		printf("\n") ; 
		}
}

/* mat1.txtのデータ */
 1   3
-1   2
 3   1
 2  -2
 0   1
 1   3
 1   0  -1
 2   1   2



/* Program 10-2 (分散共分散を計算するプログラム)　*/

#include <stdio.h>
#include <math.h>
#define NDATA 5
#define K 3
void inputdata(FILE *fdata, double x[ ][K],int m, int n);
void devcalc(double w[ ][K], double wdev[ ][K],int m, int n);
void matprod(double x[ ][NDATA], double y[ ][K], double z[ ][K], int m, int k, int n);
void mattrans(double x[ ][K], double y[ ][NDATA], int m, int n);
void out_res(double a[ ][K], int m, int n);

void main( )
{ double w[NDATA][K],wdev[NDATA][K],wdevt[K][NDATA],
wprod[K][K],s[K][K],r[K][K];
   int i,j;
FILE *fdata;
  fdata=fopen("mat2.txt","r");
inputdata(fdata,w,NDATA,K);
devcalc(w,wdev,NDATA,K);
  mattrans(wdev,wdevt,NDATA,K);
  matprod(wdevt,wdev,wprod,K,NDATA,K);
  for(i=0; i<=K-1; i++)
	  for(j=0; j<=K-1; j++)
      {s[i][j]=wprod[i][j]/(NDATA-1);
	  r[i][j]=wprod[i][j]/sqrt(wprod[i][i]*wprod[j][j]);}
	printf("  Variances and Covariances \n");
    out_res(s,K,K);
	printf("  Correlation Coefficients \n");
    out_res(r,K,K);
   fclose(fdata);
}

void inputdata(FILE *fdata,double x[ ][K],int m, int n)
{int i,j;
double d;
  for(i=0; i<= m-1; i++)
	  for(j=0; j<= n-1; j++)
	  {fscanf(fdata, "%lf", &d);
	  x[i][j]=d;
	  }
}

void devcalc(double w[ ][K], double wdev[ ][K],int m, int n)
{int i,j;
 double a;
 for(j=0; j<=n-1; j++)
 { a=0;
   for(i=0; i<=m-1; i++)a=a+w[i][j];
   a=a/m;
   for(i=0; i<=m-1; i++)wdev[i][j]=w[i][j]-a;
 }
}
 

void matprod(double x[ ][NDATA], double y[ ][K], double z[ ][K], int m, int k, int n)
{int i,j,r;
 double a;
for(i=0; i<= m-1; i++)
for(j=0; j<= n-1; j++)
{   a=0;
	for(r=0; r<=k-1; r++) a=a+x[i][r]*y[r][j];
    z[i][j]=a;
}
}

void mattrans(double x[ ][K], double y[ ][NDATA], int m, int n)
{int i,j;
for(i=0; i<=m-1; i++)
for(j=0; j<=n-1; j++)
y[j][i]=x[i][j];
}


void out_res(double x[ ][K], int m, int n)
{int i,j;
   for(i=0; i<= m-1; i++)
		{for(j=0; j<= n-1; j++)printf(" %10.4f ", x[i][j]);
		printf("\n") ; 
		}
}


/* mat2.txtのデータ */
 1  2  2
 0 -2 -1
 2  1  3
 1  3  2
-2  4  1


/* 11章のプログラム　*/

/* Program 11-1 (行列式を計算するプログラム) */

#include <stdio.h>
#include <math.h>
#define N 4
double mdeterm(double a[ ][N], int n);
void mcopy(double a[ ][N], double b[ ][N], int m, int n);
double mdet_row_trans(double b[ ][N], int n, int i, double d);

void main( )
{double a[ ][N]={{1,1,0,2},{3,-1,1,0},{0,2,2,1},{2,0,1,3}},d;
  d = mdeterm(a, N);
  printf(" |A| = %10.4f\n",d);
}

double mdeterm(double a[ ][N], int n)
{double b[N][N],d;
  int i;
  d =1;
  mcopy(a, b, n, n);
  for(i=0; i<=n-1; i++) 
  {d=mdet_row_trans(b,n,i,d);
   if(d==0)break;  
  }
  return d;
}

void mcopy(double a[ ][N], double b[][N], int m, int n)
{
	int i,j;
   for(i=0; i<=m-1;i++)
     for(j=0; j<=n-1;j++)b[i][j] = a[i][j];
}

double mdet_row_trans(double b[ ][N], int n, int i, double d)
{int i1, j, k, change;
 double b1, b2;
 i1=i-1; 
 do{
  i1++;
  b1= b[i1][i];
 }
 while(b1==0 && i1<n-1);
  
 if(b1== 0)d = 0;
 else{
 change=1;
   if(i != i1)
   {change = -1;
   for(j = i; j<=n-1; j++){
	   b2 = b[i][j];
	   b[i][j]=b[i1][j];
	   b[i1][j] = b2;}
   }
   b1 = b[i][i];
d = b1*d*change;

  for(j = i; j<= n-1; j++)b[i][j] = b[i][j]/b1;
  for( j = i+1;j<= n-1;j++){
  b1 = b[j][i];
  for(k = i; k<= n-1; k++)b[j][k] = b[j][k] - b1 * b[i][k];
 }}
  return d;
}


/* Program 11-2 (逆行列計算のプログラム) */

#include <stdio.h>
#include <math.h>
#define N 4
#define N2 2*N
double minv(double a[ ][N], double ainv[ ][N], int n);
double minv_row_trans(double b[ ][N2], int n,int n2,int i);
void out_res(double a[ ][N], int m, int n);

void main( )
{double a[ ][N]={{1,1,0,2},{3,-1,1,0},{0,2,2,1},{2,0,1,3}},
   ainv[N][N],d;
   d=minv(a,ainv,N);
   
   if(d==0)printf("Singular Matrix! \n");
   else{
   printf("Inverse Matrix \n");
    out_res(ainv,N,N);
   }
}

double minv(double a[ ][N], double ainv[ ][N], int n)
{
	double b[N][N2],d;
	int i,j,n2;
n2 = 2 * n;

  for(i = 0; i<=n-1; i++)
  {for(j = 0; j<=n-1; j++)b[i][j] = a[i][j];
    for( j = n; j<= n2-1; j++)b[i][j] = 0;
    b[i][i + n] = 1;

  }

  for(i = 0; i<= n-1; i++){
	  d=minv_row_trans(b, n, n2, i);
	  if(d==0) break;}

    for(i = 0; i<=n-1; i++)
    for(j = 0; j<=n-1; j++)ainv[i][j] = b[i][j + n];
 return d;
}

double minv_row_trans(double b[ ][N2], int n,int n2,int i)
{
int i1,j,k;
double b1,b2;
   i1 = i - 1;
   do
   {i1++;
    b1 = b[i1][i];
	
   }while(b1== 0 && i1 < n-1);
     
   if(b1!=0){
      if(i1 != i)
	  for(j = 0; j<= n2-1; j++)
	  {b2 = b[i][j];
      b[i][j] = b[i1][j];
      b[i1][j] = b2;}

	  b1 = b[i][i];
      for(j = 0; j<= n2-1; j++)b[i][j] = b[i][j] / b1;

    for(j = 0; j<= n-1; j++){
	  if(i !=j)
	  {b1 = b[j][i];
      for (k = 0; k<= n2-1; k++)b[j][k] = b[j][k] - b1 * b[i][k];
	  }
	}
	}
return b1;}

void out_res(double x[ ][N], int m, int n)
{int i,j;
   for(i=0; i<= m-1; i++)
		{for(j=0; j<= n-1; j++)printf(" %12.5f ", x[i][j]);
		printf("\n") ; 
		}
}



/* 12章のプログラム　*/

/* Program 12-1 (最小二乗推定量を計算するプログラム)　*/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>

#define K 3
#define N 6
#define K2 2*K
#define PI 3.141592654

void inputdata(FILE *fdata, double x[ ][K], double y[ ], int m, int n);
void matprod(double *x, double *y, double *xy, int m, int r,int n);
void mattrans(double *x, double *y, int m, int n);
void out_res(FILE *fout, double *x, int m, int n);
double ls_calc(double x[ ][K],double y[ ], double beta[ ][3],double yfit[],
		 double e[ ],double *ssr,double *s2,double var[ ][K], int m, int n);
double minv(double a[ ][K], double ainv[ ][K], int n);
double minv_row_trans(double b[ ][K2], int n,int n2,int i);
void ystat(double y[ ], double *ymean, double *ystd, double *ydevsq, int n);
void reg_stat(double e[ ],double ssr,double *r2, double *adjr2, double ydevsq,
		  double *logL, double *aic, double *f, double *dw, int n, int k);
double dw_calc(double e[ ], int n, double ssr);


void main( )
{ double x[N][K], y[N], beta[K][3], yfit[N],e[N],var[K][K], 
         ssr,s2,ymean,ystd,ydevsq,r2,adjr2,logL,aic,f,dw,d;
FILE *fdata,*fout;
  if((fdata=fopen("reg1.txt","r"))==NULL)
  {printf("File does not exist!\n");
  exit(1);}

  fout=fopen("res1.out","w");
  inputdata(fdata,x,y,N,K);
  d=ls_calc(x,y,beta,yfit,e, &ssr,&s2,var,N,K);
  
  if(d != 0){
   ystat(y,&ymean, &ystd, &ydevsq, N);
   reg_stat(e,ssr,&r2,&adjr2,ydevsq,&logL,&aic,&f,&dw,N,K);
   
    fprintf(fout,"Estimation Results: Least Squres Method \n\n");
    fprintf(fout,"Number of Observations:%5d     Degree of Freedom:%5d\n",N,N-K);
    fprintf(fout,"Mean of Dependent Varible:      %12.5f \n", ymean); 
    fprintf(fout,"Std. Dev. of Dependent Variable:%12.5f \n",ystd); 
    fprintf(fout,"Sum of Squared Residuals:       %12.5f\n",ssr);
    fprintf(fout,"Std. Error of Regression:       %12.5f\n", sqrt(s2));
    fprintf(fout,"R2:  %12.5f    Adusted R2:%12.5f \n", r2,adjr2);
    fprintf(fout,"LogL:%12.5f    AIC:       %12.5f \n", logL, aic);
    fprintf(fout,"F-Statistic:%12.5f  Durbin-Watson Statistic:%12.5f\n\n",f,dw);  
    
	 fprintf(fout,"  Coefficient    Std. Error   t-statistic \n");
     out_res(fout,beta[0],K,3);
     fprintf(fout,"\n  Variance-Covariance Matrix \n");
     out_res(fout,var[0],K,K);
  }
  else printf("Cannot be estimated! Singular Matrix! \n");
 
   fclose(fdata);
   fclose(fout);
}

void inputdata(FILE *fdata,double x[ ][K],double y[ ],int m, int n)
{int i,j;
double d;
  for(i=0; i<= m-1; i++)
	   {fscanf(fdata, "%lf", &d);
        y[i]=d;
		x[i][0]=1;
	  for(j=1; j<= n-1; j++)
	  {fscanf(fdata, "%lf", &d);
	  x[i][j]=d;
	  }
  }
}

double ls_calc(double x[ ][K],double y[ ], double beta[ ][3],double yfit[],
	  double e[ ],double *ssr,double *s2,double var[ ][K], int n, int k)
{int i,j;
 double xt[K][N], xtx[K][K], xtxinv[K][K], xty[K],beta1[K],d;
  mattrans(x[0],xt[0],n,k);
  matprod(xt[0],x[0],xtx[0],k,n,k);
  d=minv(xtx,xtxinv,k);
  
 if(d != 0)
  {
  matprod(xt[0],y,xty,k,n,1);
  matprod(xtxinv[0],xty,beta1,k,k,1);
    matprod(x[0],beta1,yfit,n,k,1);
	for(i=0; i<=n-1; i++)e[i]=y[i]-yfit[i];
	matprod(e,e,ssr,1,n,1);
	   *s2=*ssr/(n-k);
    for(i=0; i<=k-1; i++)
	   {beta[i][0]=beta1[i];
	    beta[i][1]=sqrt(*s2*xtxinv[i][i]);
		beta[i][2]=beta[i][0]/beta[i][1];
	    for(j=0; j<=k-1; j++) var[i][j]=*s2*xtxinv[i][j];
	}
  }
 return d;
}

void matprod(double *x,double *y,double *xy, int n, int r, int m)
{int i1,i2,i3;
  double a,*x1,*y1;
  for(i1=0; i1<=n-1; i1++){
  for(i2=0; i2<=m-1; i2++)
  {a=0;
   x1=x+i1*r;
   y1=y+i2;
    for(i3=0; i3<=r-1; i3++){
	  a=a+(*x1)*(*y1);
	   x1++;
      y1=y1+m;}
  *xy=a;
  xy++;}     
  }
}


void mattrans(double *x, double *y, int m, int n)
{int i,j,k1,k2;
 double *x1, *y1;
   for(i=0; i<=m-1; i++)
     for(j=0; j<=n-1; j++){
     k1=i*n + j;
     k2=j*m + i;
     x1=x+k1;
     y1=y+k2;
     *y1=*x1;
	 }
}

double minv(double a[ ][K], double ainv[ ][K], int n)
{ double b[K][K2],d;
  int i,j,n2;
  n2 = 2 * n;
  for(i = 0; i<=n-1; i++)
  {for(j = 0; j<=n-1; j++)b[i][j] = a[i][j];
  for( j = n; j<= n2-1; j++)b[i][j] = 0;
  b[i][i + n] = 1;
  }
  for(i = 0; i<= n-1; i++){
	  d=minv_row_trans(b, n, n2, i);
	  if(d==0) break;}
 for(i = 0; i<=n-1; i++)
	 for(j = 0; j<=n-1; j++)ainv[i][j] = b[i][j + n];
	 return d;
}

double minv_row_trans(double b[ ][K2], int n,int n2,int i)
{int i1,j,k;double b1,b2;i1 = i - 1;
   do
   {i1++;
    b1 = b[i1][i];
   }while(b1== 0 && i1 < n-1);
	if(b1!=0)
	{if(i1 != i)
	for(j = 0; j<= n2-1; j++)
	{b2 = b[i][j];
	b[i][j] = b[i1][j];
	b[i1][j] = b2;}
	b1 = b[i][i];
	for(j = 0; j<= n2-1; j++)b[i][j] = b[i][j] / b1;
       for(j = 0; j<= n-1; j++)
	   {if(i!=j){b1 = b[j][i];
	   for (k = 0; k<= n2-1; k++)
		b[j][k] = b[j][k] - b1 * b[i][k];}
	}
}
return b1;}

void ystat(double y[ ], double *ymean, double *ystd, double *ydevsq, int n)
{int i;
  double a; 
   a=0;
  for(i=0; i<=n-1; i++) a=a+y[i];
  *ymean=a/n;
   a=0;
   for(i=0; i<=n-1; i++) a=a+pow((y[i]-*ymean),2);
   *ydevsq=a;
   *ystd=sqrt(a/(n-1));
}

void reg_stat(double e[ ],double ssr,double *r2, double *adjr2, double ydevsq,
		double *logL, double *aic, double *f, double *dw, int n, int k)
{
	*r2=1-ssr/ydevsq;
	*adjr2=1-(ssr/(n-k))/(ydevsq/(n-1));
	*logL=-(n/2)*(1+log(2*PI)+log(ssr/n));
	*aic=-2*(*logL)+2*k;
	*f=((ydevsq-ssr)/(k-1))/(ssr/(n-k));
	*dw=dw_calc(e,n,ssr);
}

double dw_calc(double e[ ], int n, double ssr)
{double a;
 int i;
 a=0;
 for(i=0; i<=n-2; i++)a=a+pow((e[i+1]-e[i]),2);
 return a/ssr;
}


void out_res(FILE *fout,double *x, int m, int n)
{int i,j;
for(i=0; i<= m-1; i++)
	{
	for(j=0; j<= n-1; j++)
	{fprintf(fout," %12.6f ", *x);
	x++;}
	fprintf(fout,"\n");
}
} 

/*  reg1.txtのデータ　*/

2  2  4
3  2  3
1  3  3
4  4  1
5  4  2
4  5  1


/*　13章のプログラム　*/

/* Program 13-1 (ウィルコクスンの順位和検定を行うプログラム) */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define N1 5
#define N2 6
#define N N1+N2
#define K 3
void inputdata(FILE *fdata, double z[ ][K], int n1, int n2);
void sort2(double z[ ][K], int n, int k, int m, int order);
void same_val(double z[ ][K], int n);
double norminv(double u);
void wrank_sum_percent(int n1, int n2, double per[ ][3]);

void main( )
{ double z[N][K],rsum,per[8][3];
  int i,j;
  FILE *fdata;
    per[0][0] = 0.005;
    per[1][0] = 0.01;
    per[2][0] = 0.025;
    per[3][0] = 0.05;
    per[4][0] = 0.1;
    per[5][0] = 0.2;
    per[6][0] = 0.3;
    per[7][0] = 0.4;
  if((fdata=fopen("wrank1.txt","r"))==NULL)
  {printf("File does not exist!\n");
  exit(1);}
  inputdata(fdata,z,N1,N2);
     sort2(z, N, K, 1,0);
	 for(i=0; i<=N-1; i++)z[i][2]=i+1;
	 same_val(z, N);
     rsum=0;
     for(i=0; i<=N-1; i++) if(z[i][1]==1)rsum=rsum+z[i][2];
	 wrank_sum_percent(N1, N2, per);

   printf("Wilcoxson Test\n");
   printf("n  = %5d\nn1 = %5d \nn2 = %5d\n", N, N1,N2);
   printf("Sum of Ranks = %12.3f\n\n",rsum);
   printf("Percentile: \n");
   printf("       alpha       Lower       Upper\n");
   for(i=0; i<=7; i++)
		{for(j=0; j<=2; j++)printf("%12.3f", per[i][j]);
		 printf("\n");}
   fclose(fdata);
}

void inputdata(FILE *fdata,double z[ ][K],int n1, int n2)
{int i,n;
double d1,d2;
n=n1+n2;
if(n2 >= n1){ 
  for(i=0; i<= 2*(n1-1); i=i+2)
	   {fscanf(fdata, "%lf %lf", &d1, &d2);
        z[i][0]=d1;
		z[i][1]=1;
		z[i+1][0]=d2;
		z[i+1][1]=2;
		}
  for(i=2*n1; i<=n-1; i++)
	  {fscanf(fdata, "%lf", &d1);
	  z[i][0]=d1;
	 z[i][1]=2;}
}
if(n1>n2){ 
  for(i=0; i<= 2*(n2-1); i=i+2)
	   {fscanf(fdata, "%lf %lf", &d1, &d2);
        z[i][0]=d1;
		z[i][1]=2;
		z[i+1][0]=d2;
		z[i+1][1]=1;}
  for(i=2*n2; i<=n-1; i++)
	  {fscanf(fdata, "%lf", &d1);
	  z[i][0]=d1;
	  z[i][1]=2;}
}
 }


void sort2(double z[ ][K], int n, int k, int m, int order)
{
int change, i,j,n1;
double a;
m=m-1;
do
{
change = 0;
  for(i=0; i<=n - 2; i++)
  {
     if(z[i][m] > z[i + 1][m])
	 {change = 1;
		for(j = 0; j<= k-1; j++)
		{a = z[i][j];
		z[i][j] = z[i + 1][j];
		z[i+1][j] = a;}
	 }
  }
}while(change ==1);
if(order == 1)
{n1 = n/2;
	for(i = 0; i <= n1-1; i++)
	{for(j = 0; j <= k-1; j++)
		{a = z[i][j];
		z[i][j] = z[n-1-i][j];
		z[n-1-i][j] = a;}
	}
}
}

void same_val(double z[ ][K], int n)
{
int same_ord[N], ord_start[N][2],i, i1, i2, j, j1, j2,sd1;
double a;
  for(i = 0; i <= n - 2; i++)
  {if(z[i][0] == z[i + 1][0])same_ord[i] = 1;
  else same_ord[i] = 0;
  }    
same_ord[n-1] = 0;
j1 = 0;
for(i = 0; i<= n - 1; i++)
{if(i==0)sd1 =0;
else sd1= same_ord[i - 1];
  if(sd1== 0 && same_ord[i] == 1)
	{ord_start[j1][0] = i;
	 j2=0;
		do{
		j2++;
		}while(same_ord[i+j2] ==1);
	ord_start[j1][1] = j2;
	j1++;
  }
}
if(j1 >= 1) 
	{for(i = 0; i<= j1-1; i++)
		{a = 0;
		i1 = ord_start[i][0];
		i2 = ord_start[i][1];
		for(j = i1; j <= i1+i2; j++)a = a + z[j][2];
        for(j = i1; j <= i1+i2; j++)z[j][2] = a / (i2 + 1);
	}
}
}

double norminv(double u)
{double a,b,w,z;
  z=-log(4*u*(1-u));
  a=z*(2.0611786-5.7262204/(z+11.640595));
  w=sqrt(a);
  if(u<0.5)b=-w;
  else b=w;
return b;
}

void wrank_sum_percent(int n1, int n2, double per[ ][3])
{ double mean, var, s, p, a;
  int n, i;
  n = n1 + n2;
  mean = 1.0*n1 *(n + 1) / 2.0;
  var = 1.0*n1 * n2 * (n + 1) / 12.0;
  s = sqrt(var);
    for(i = 0; i <= 7; i++)
	{p = per[i][0];
     a = -norminv(p) * s;
     per[i][1] = mean - a;
     per[i][2] = mean + a;
	}
}


/* wrank1.txtのデータ */

 3.2  3.2
 3.4  4.5
 5.6  2.5
 4.1  6.2
 2.9  7.1
      4.3


/* Program 13-2 (ウィルコクスンの符号付順位検定)　*/

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

#define N 10
#define K 4
void inputdata(FILE *fdata, double z[ ][K], int n);
void sort2(double z[ ][K], int n, int k, int m, int order);
void same_val2(double z[ ][K], int n);
double norminv(double u);
void wrank_sign_percent(int n1, double per[ ][3]);

void main( )
{ double z[N][K],rsum,per[8][3];
  int i,j,n1;
  FILE *fdata;
    per[0][0] = 0.005;
    per[1][0] = 0.01;
    per[2][0] = 0.025;
    per[3][0] = 0.05;
    per[4][0] = 0.1;
    per[5][0] = 0.2;
    per[6][0] = 0.3;
    per[7][0] = 0.4;
  if((fdata=fopen("wrsum1.txt","r"))==NULL)
  {printf("File does not exist!\n");
  exit(1);}
  inputdata(fdata,z,N);
	sort2(z, N, K, 3,0);
	for(i=0; i<=N-1; i++)z[i][3]=i+1;
	same_val2(z, N);
rsum=0;
	n1=0;
	for(i=0; i<=N-1; i++) 
		{if(z[i][0] > z[i][1])rsum=rsum+z[i][3];
		if(z[i][0] != z[i][1])n1++;
		}
   printf("Wilcoxson Signed Rank Test\n");
   printf("n  = %5d  \nn1 = %5d \n", N, n1);
   printf("Sum of Ranks = %12.3f\n\n",rsum);
   if(n1>20){
   wrank_sign_percent(n1, per);
   printf("Percentile: \n");
   printf("       alpha       Lower       Upper\n");
   for(i=0; i<=7; i++)
		{for(j=0; j<=2; j++)printf("%12.3f", per[i][j]);
		 printf("\n");}
   }
   fclose(fdata);
}

void inputdata(FILE *fdata,double z[ ][K],int n)
{int i;
double d1,d2;
  for(i=0; i<= n-1; i++)
	   {fscanf(fdata, "%lf %lf", &d1, &d2);
        z[i][0]=d1;
		z[i][1]=d2;
		z[i][2]=fabs(d1-d2);
		if(z[i][2]==0)z[i][2]=1.111e100;
		}
}


void sort2(double z[ ][K], int n, int k, int m, int order)
{
int change, i,j,n1;
double a;
m=m-1;
do
{
change = 0;
  for(i=0; i<=n - 2; i++)
  {
     if(z[i][m] > z[i + 1][m])
	 {change = 1;
		for(j = 0; j<= k-1; j++)
		{a = z[i][j];
		z[i][j] = z[i + 1][j];
		z[i+1][j] = a;}
	 }
  }
}while(change ==1);
if(order == 1)
{n1 = n/2;
	for(i = 0; i <= n1-1; i++)
	{for(j = 0; j <= k-1; j++)
		{a = z[i][j];
		z[i][j] = z[n-1-i][j];
		z[n-1-i][j] = a;}
	}
}
}

void same_val2(double z[ ][K], int n)
{
int same_ord[N], ord_start[N][2],i,i1,i2,j,j1,j2,sd;
double a;
  for(i = 0; i <= n-2; i++)
  {if((z[i][2] == z[i + 1][2]) && (z[i][0] != z[i][1]))
	same_ord[i] = 1;
    else same_ord[i] = 0;
  }

same_ord[n-1] = 0;
j1=0;
  for(i = 0; i <= n - 1; i++)
  {if(i==0)sd=0;
	else sd=same_ord[i-1];
      if((sd ==0) && (same_ord[i] == 1))
		{ord_start[j1][0] = i;
			j2 = 0;
			do{
			j2++;
			}while(same_ord[i + j2] ==1);
       ord_start[j1][1]=j2;
	   j1++;
	  }
  }
if(j1 >= 1)
{
  for(i = 0; i<= j1-1; i++)
  {a = 0;
  i1 = ord_start[i][0];
  i2 = ord_start[i][1];
    for(j = i1; j<= i1 + i2; j++)a = a + z[j][3];
    for(j = i1; j<= i1 + i2; j++)z[j][3] = a / (i2 + 1);
  }
}
}

double norminv(double u)
{double a,b,w,z;
  z=-log(4*u*(1-u));
  a=z*(2.0611786-5.7262204/(z+11.640595));
  w=sqrt(a);
  if(u<0.5)b=-w;
  else b=w;
return b;
}

void wrank_sign_percent(int n,double per[ ][3])
{ double mean, var, s, p, a;
  int i;
  mean = n * (n + 1) / 4;
  var = n * (n + 1) * (2 * n + 1) / 24;
  s = sqrt(var);
    for(i = 0; i <= 7; i++)
	{p = per[i][0];
     a = -norminv(p) * s;
     per[i][1] = mean - a;
     per[i][2] = mean + a;
	}
}

/* wrsum1.txtのデータ　*/

 3.2  2.1
 3.4  4.5
 5.6  2.5
 7.2  6.2
 2.9  7.1
 7.0  4.3
 4.2  4.2
 3.6  4.2
 6.1  2.5
 5.3  7.3