数量化Ｉ類プログラムquant1.c
/* 多変量統計解析法 第4章 §２数量化Ⅰ類（pp.139-151)
(河口至商、多変量解析入門Ⅰ、森北出版、1973、pp.95-107)
quant1.c 数量化Ⅰ類, 1994.01.20, 1998.1.29, 1999.2.18 z.sekiya */
/* A. プログラム仕様 program spec (quant1.c)
1) 入力 input:
title：タイトル, n:サンプル数sample size, ns:アイテム数items ,
kn[]:カテゴり数 categorys
y,jz[n][ns]
2) 処理 process:
x：データ行列, v(<= x'x)：生の積和行列 , 掃出しsweep out,
カテゴり毎の標準化normalize category weight
x: データ行列（ダミー変数dummy variable,外的基準y,定数1）
3) 出力 output:
外的基準y,jz, データ行列x, 積和行列v, v^-1, カテゴリスコア, 範囲
予測値、誤差、重相関係数、偏相関係数
B. 処理手順 plan
1) データファイル名の入力 Input file name
2) データの読み込み read from file ( title, parameters, data )
3) データ行列ｘへのセット set x
4) 積和行列、分散の計算 calculate v, print v
5) 掃き出し sweep s, print s
6) カテゴリ毎の標準化 normalise category weight, print b, bh
7) 予測値と誤差の計算・表示 y, yhat, e
8) 標本行列 X、生の積和行列 rs(x'x),積和行列 s(x'x)の計算
9) 相関行列 R(x'x), 偏相関係数 sweep out for pr
1. (compile and link) >cl quant1.c sweep.obj
2. (Execution) >quant1
(file name needed for data input ?) file_name
*/
#include
#include
#include
#define LA 12
#define LN 30
void sweep(double a[][LA], int k, int ii); // 掃き出し
void correl(double s[][LA], int jj, double r[][LA]); // 相関係数の計算
void matprt(char name[], char title[], double a[][LA], int im, int jm);// 行列表示
matprod(double x[][LA],int n, int nnp2, double v[][LA]) // 生の積和行列の計算
{
int i, j, l; double vjl;
for(j=1; j<=nnp2; j++) {
for(l=1; l<=nnp2; l++) {
vjl=0.0;
for(i=1;i<=n; i++)
vjl=vjl+x[i][j]*x[i][l];
v[j][l]=vjl;
}
}
}
main() // 数量化Ⅰ類,main関数
{
char fname[30], title[60], *pa, aster[5];
int i, j, k,kk, l, m, n, ns, kn[LA], kc, nn, nnp1, is, jz[LN][LA],
nnp2, jcz[LA][LA], freq[LA], freq1[LA], nsp1, nsp2;
double yb, y[LN], yy[LN],e[LN], x[LN][LA], v[LA][LA], yhat,
vjl, bh[LA], xb[LA], aij[LA][LA], an[LA][LA], sumc[LA],
anmax[LA],anmin[LA],anrange,y2,yh2, s[LA][LA],
pr[LA], r[LA][LA];
FILE *fp1;
long nowtime;
time(&nowtime);
pa = ctime(&nowtime);
*(pa+24)='\0';
printf("<%s>\n", pa);
/* 1) データファイル名の入力 Input file name */
do {
printf("<< 数量化Ⅰ類(quant1) k134084 関谷 >> Input ファイル名:");scanf("%s", fname);
fp1 = fopen(fname,"r");
} while(fp1 ==NULL) ;
/* 2) データの読み込み read from file (title, parameters, data ) */
/* 3) データ行列ｘへのセット set x */
fscanf(fp1, "%s %d %d",title, &n, &ns);
for(i=1; i<=ns; i++) fscanf(fp1, "%d", &kn[i]);
printf("\nタイトル<%s>\nサンプル数(n):%d アイテム数(ns):%d \n", title, n, ns);
printf(" アイテム毎のカテゴリ数(kn[i]) ");
kc=0;
for(i=1; i<=ns; i++) {
printf(" kn[%d]:%d ", i, kn[i]);
kc += kn[i];
}
nn = kc - ns + 1; nnp1 = nn + 1; nnp2=nnp1+1;
printf("\n ダミー変数の数(kc += kn[i]):%d\n", kc);
printf("\n 方程式の実際の変数の数(nn=kc-ns+1):%d\n", nn);
printf(" 変数添え字への変換テーブルjcz[アイテム][カテゴリ]\n");
k=0;
for(i=1;i<=ns;i++)
for(j=1;j<=kn[i]; j++) { // これでデータ行列Ｘをセットする
if(j==1 && i>1) jcz[i][j] = 0; else jcz[i][j] = ++k;
printf(" jcz[%d][%d]:%d ",i,j, jcz[i][j]);
}
printf("\n 方程式の実際の変数の数の確認(k):%d\n", k);
printf("\nサンプルデータの読み込み\n");
printf(" サンプル番号(i) 目的変数(y[i]) アイテムjのカテゴリ番号jz[i][1] ");
for (j=2; j<=ns; j++) printf(" jz[i][%d] ", j);
for(i=1; i<=n; i++){
for(k=1;k<=nn; k++)
x[i][k]=0.0;
x[i][nnp2]=1.0;
fscanf(fp1, "%d %lf", &kk, &x[i][nnp1]); y[i]= x[i][nnp1];
printf("\n i:%2d y[i]: %f ", kk, x[i][nnp1]);
for(j=1;j<=ns;j++){
fscanf(fp1, "%d", &jz[i][j]);printf(" jz[%2d][%d]:%d ", i, j, jz[i][j]);
x[i][(jcz[j][(jz[i][j])])]=1.0; // データ行列へのセット
}
}
printf("\n\nダミー変数への代入 x[i][(jcz[j][(jz[i][j])])]=1.0;\n");
printf("アイテムNo "); // 行の見出し（データ行列 X）
for (i=1;i<=ns; i++){
for (j=1; j<=kn[i]; j++){
if (i>1 && j ==1) continue;
printf(" %d ", i);
}
}
printf(" 目的変数y ダミー変数\nカテゴリNo ");
for (i=1;i<=ns; i++){
for (j=1; j<=kn[i]; j++){
if (i>1 && j ==1) continue;
printf(" %d ", j);
}
}
matprt("(quant1.c)データ行列 X seted", title, x, n, nnp2);
/* 4) 積和行列、分散の計算 v=x'x */
matprod(x,n,nnp2, v);
matprt("積和行列 v(x'x)", title, v, nnp2, nnp2);
for(i=1;i<=nn;i++) {
freq[i]=v[nnp2][i]; // カテゴリの度数
xb[i]=freq[i]/(float)n; // カテゴリの相対度数
}
yb=v[nnp2][nnp1]/n; // 目的関数の平均,定数項である
/* 5) 掃き出し sweep out */
for (k=1; k<=nn; k++) {
sweep(v, k, nnp2);
}
matprt("掃き出し後の v-1 ", title, v, nnp2, nnp2);
/* 6) カテゴリ毎の標準化 normalise category weight, print b, bh */
k=1;
for(i=1;i<=ns; i++){
sumc[i]=0.0; freq1[i]=n; anmax[i]=-10000.0; anmin[i]=10000.0;
for(j=1;j<=kn[i];j++) {
if(i>1 && j==1) {
an[i][j]=0.0; // アイテムiカテゴリｊへの係数
} else {
an[i][j]=v[k++][nnp1]; // 正規方程式の解
sumc[i] += an[i][j]*xb[k-1]; // 相対度数との積和（アイテム毎）
freq1[i] -= freq[k-1]; // アイテム毎のカテゴリ１の度数
}
aij[i][j]=an[i][j]; // 偏相関係数のため
if(an[i][j] >anmax[i]) anmax[i]=an[i][j]; // アイテム内最大の係数
if(an[i][j] 1) {
printf(" %d %2d %10.2f\n", j, freq1[i], an[i][j]);
}else
printf(" %d %2d %10.2f\n", j, freq[k++], an[i][j]);
}
}
printf(" 定数項（ｙの平均:yb）: %8.2f\n", yb);
/* 7) 予測値と誤差の計算・表示 y, yhat, e、print yhat, e */
printf("\n<< 予測誤差など[ %s ] >>>\n No. 観測値 予測値 誤差 \n", title);
y2=yh2=0.0;
for(i=1;i<=n;i++){
yhat = yb;
for(j=1; j<=ns; j++) {
k=jz[i][j];
yhat = yhat + an[j][k];
}
e[i]= y[i]-yhat;
printf(" %2d %6.2f %6.2f %6.2f ", i, y[i], yhat, e[i]);
y2 += (y[i]-yb)*(y[i]-yb);
yh2 += (yhat-yb)*(yhat-yb);
for(j=1;j<=ns;j++)
printf(" %2d", jz[i][j]);
printf("\n");
}
printf("\n 重相関係数 Ry12...m: %6.3f\n", sqrt(yh2/y2));
/* 8) 標本行列 X、生の積和行列 rs(x'x),積和行列 s(x'x)の計算 */
nsp1=ns+1; nsp2=ns+2;
for(i=1; i<=n; i++){
for(j=1; j<=ns; j++)
x[i][j] = aij[j][jz[i][j]]; // アイテム得点
x[i][nsp1] = y[i]; // 目的関数
x[i][nsp2] = 1.0; // 定数項
}
matprt("標本行列 X(アイテム得点1,2,..,m,y,1.0)",title, x, n, nsp2);
matprod(x,n,nsp2, s);
matprt(" 生の積和行列 rs(x'x)", title, s, nsp2, nsp2);
sweep(s,nsp2,nsp2);
matprt(" 積和行列 s(x'x)", title, s, nsp2, nsp2);
/* 9) 相関行列 R(x'x), 偏相関係数 sweep out for pr */
correl(s,nsp1,r);
matprt(" 相関行列 R(x'x)", title, r, nsp1, nsp1);
for (k=1; k<=nsp1; k++)
sweep(r, k, nsp1);
matprt("掃き出し後 R-1 ", title, r, nsp1, nsp1);
for(i=1;i<=ns;i++)
pr[i]=-r[i][nsp1]/sqrt(r[i][i]*r[nsp1][nsp1]);
printf("\n<偏相関係数> \n");
for(i=1; i<=ns; i++)
printf(" pr(y,%d):%7.3f\n", i, pr[i]);
printf("\n\n");
}
/* 以下は、実行例(qut1.out)を参照のこと
*/