| 情報学部 | 菅沼ホーム | 目次 | 索引 |
/****************************/
/* 正準相関分析 */
/* coded by Y.Suganuma */
/****************************/
#include <stdio.h>
#include <math.h>
int gauss(double **, int, int, double);
int canonical(int, int, int, double **, double *, double **, double, double *, int, int);
int power(int, int, int, int, double, double **, double **,
double **, double *, double **, double *, double *, double *);
int main()
{
double **x, *ryz, **ab, *v0;
int i1, i2, q, r, s, n, sw;
scanf("%d %d %d", &r, &s, &n); // 各組の変数の数とデータの数
q = r + s;
ryz = new double [r];
v0 = new double [r];
x = new double * [q];
ab = new double * [r];
for (i1 = 0; i1 < q; i1++)
x[i1] = new double [n];
for (i1 = 0; i1 < r; i1++) {
ab[i1] = new double [q];
v0[i1] = 1.0;
}
for (i1 = 0; i1 < n; i1++) { // データ
for (i2 = 0; i2 < q; i2++)
scanf("%lf", &x[i2][i1]);
}
sw = canonical(r, s, n, x, ryz, ab, 1.0e-10, v0, 15, 200);
if (sw > 0) {
for (i1 = 0; i1 < sw; i1++) {
printf("相関係数 %f\n", ryz[i1]);
printf(" a");
for (i2 = 0; i2 < r; i2++)
printf(" %f", ab[i1][i2]);
printf("\n");
printf(" b");
for (i2 = 0; i2 < s; i2++)
printf(" %f", ab[i1][r+i2]);
printf("\n");
}
}
else
printf("***error 解を求めることができませんでした\n");
for (i1 = 0; i1 < q; i1++)
delete [] x[i1];
for (i1 = 0; i1 < r; i1++)
delete [] ab[i1];
delete [] x;
delete [] ab;
delete [] ryz;
delete [] v0;
return 0;
}
/***********************************/
/* 正準相関分析 */
/* r,s : 各組における変数の数 */
/* n : データの数 */
/* x : データ */
/* ryz : 相関係数 */
/* ab : 各組の係数（a,bの順） */
/* eps : 正則性を判定する規準 */
/* v0 : 固有ベクトルの初期値 */
/* m : 丸め誤差の修正回数 */
/* ct : 最大繰り返し回数 */
/* return : >0 : 相関係数の数 */
/* =0 : エラー */
/***********************************/
int canonical(int r, int s, int n, double **x, double *ryz, double **ab, double eps,
double *v0, int m, int ct)
{
double **A, **B, **C, **C11, **C11i, **C12, **C21, **C22, **C22i,
*mean, vv, **w, *w1, len, *v1, *v2;
int i1, i2, i3, q, sw = 0, sw1;
// 領域の確保
q = r + s;
mean = new double [q];
w1 = new double [s];
v1 = new double [r];
v2 = new double [r];
A = new double * [s];
B = new double * [r];
C = new double * [q];
C11 = new double * [r];
C11i = new double * [r];
C12 = new double * [r];
C21 = new double * [s];
C22 = new double * [s];
C22i = new double * [s];
w = new double * [s];
for (i1 = 0; i1 < q; i1++)
C[i1] = new double [q];
for (i1 = 0; i1 < r; i1++) {
B[i1] = new double [r];
C11[i1] = new double [r];
C11i[i1] = new double [r];
C12[i1] = new double [s];
}
for (i1 = 0; i1 < s; i1++) {
A[i1] = new double [s];
C21[i1] = new double [r];
C22[i1] = new double [s];
C22i[i1] = new double [s];
w[i1] = new double [2*s];
}
// 平均値の計算
for (i1 = 0; i1 < q; i1++) {
mean[i1] = 0.0;
for (i2 = 0; i2 < n; i2++)
mean[i1] += x[i1][i2];
mean[i1] /= n;
}
// 分散強分散行列の計算
for (i1 = 0; i1 < q; i1++) {
for (i2 = i1; i2 < q; i2++) {
vv = 0.0;
for (i3 = 0; i3 < n; i3++)
vv += (x[i1][i3] - mean[i1]) * (x[i2][i3] - mean[i2]);
vv /= (n - 1);
C[i1][i2] = vv;
if (i1 != i2)
C[i2][i1] = vv;
}
}
// C11, C12, C21, C22 の設定
// C12
for (i1 = 0; i1 < r; i1++) {
for (i2 = 0; i2 < s; i2++)
C12[i1][i2] = C[i1][i2+r];
}
// C21
for (i1 = 0; i1 < s; i1++) {
for (i2 = 0; i2 < r; i2++)
C21[i1][i2] = C[i1+r][i2];
}
// C11とその逆行列
for (i1 = 0; i1 < r; i1++) {
for (i2 = 0; i2 < r; i2++) {
w[i1][i2] = C[i1][i2];
C11[i1][i2] = C[i1][i2];
}
for (i2 = r; i2 < 2*r; i2++)
w[i1][i2] = 0.0;
w[i1][i1+r] = 1.0;
}
sw1 = gauss(w, r, r, 1.0e-10);
if (sw1 == 0) {
for (i1 = 0; i1 < r; i1++) {
for (i2 = 0; i2 < r; i2++)
C11i[i1][i2] = w[i1][i2+r];
}
}
else
sw = 1;
// C22とその逆行列
for (i1 = 0; i1 < s; i1++) {
for (i2 = 0; i2 < s; i2++) {
w[i1][i2] = C[i1+r][i2+r];
C22[i1][i2] = C[i1+r][i2+r];
}
for (i2 = s; i2 < 2*s; i2++)
w[i1][i2] = 0.0;
w[i1][i1+s] = 1.0;
}
sw1 = gauss(w, s, s, eps);
if (sw1 == 0) {
for (i1 = 0; i1 < s; i1++) {
for (i2 = 0; i2 < s; i2++)
C22i[i1][i2] = w[i1][i2+s];
}
}
else
sw = 1;
// 固有値λ及び固有ベクトルaの計算
if (sw == 0) {
// 行列の計算
for (i1 = 0; i1 < s; i1++) {
for (i2 = 0; i2 < r; i2++) {
A[i1][i2] = 0.0;
for (i3 = 0; i3 < s; i3++)
A[i1][i2] += C22i[i1][i3] * C21[i3][i2];
}
}
for (i1 = 0; i1 < r; i1++) {
for (i2 = 0; i2 < s; i2++) {
w[i1][i2] = 0.0;
for (i3 = 0; i3 < s; i3++)
w[i1][i2] += C12[i1][i3] * A[i3][i2];
}
}
for (i1 = 0; i1 < r; i1++) {
for (i2 = 0; i2 < r; i2++) {
A[i1][i2] = 0.0;
for (i3 = 0; i3 < r; i3++)
A[i1][i2] += C11i[i1][i3] * w[i3][i2];
}
}
// 固有値と固有ベクトル（べき乗法）
for (i1 = 0; i1 < r; i1++) {
for (i2 = 0; i2 < r; i2++)
B[i1][i2] = 0.0;
B[i1][i1] = 1.0;
}
sw = power(0, r, m, ct, eps, A, B, C, ryz, ab, v0, v1, v2);
if (sw > 0) {
for (i1 = 0; i1 < r; i1++) {
// 相関係数
ryz[i1] = sqrt(ryz[i1]);
// 大きさの調整(a)
for (i2 = 0; i2 < r; i2++) {
w1[i2] = 0.0;
for (i3 = 0; i3 < r; i3++)
w1[i2] += C11[i2][i3] * ab[i1][i3];
}
len = 0.0;
for (i2 = 0; i2 < r; i2++)
len += ab[i1][i2] * w1[i2];
len = sqrt(len);
for (i2 = 0; i2 < r; i2++)
ab[i1][i2] /= len;
// bの計算
for (i2 = 0; i2 < s; i2++) {
w1[i2] = 0.0;
for (i3 = 0; i3 < r; i3++)
w1[i2] += C21[i2][i3] * ab[i1][i3];
}
for (i2 = 0; i2 < s; i2++) {
ab[i1][i2+r] = 0.0;
for (i3 = 0; i3 < s; i3++)
ab[i1][i2+r] += C22i[i2][i3] * w1[i3];
}
for (i2 = 0; i2 < s; i2++)
ab[i1][i2+r] /= ryz[i1];
// 大きさの調整(b)
for (i2 = 0; i2 < s; i2++) {
w1[i2] = 0.0;
for (i3 = 0; i3 < s; i3++)
w1[i2] += C22[i2][i3] * ab[i1][i3+r];
}
len = 0.0;
for (i2 = 0; i2 < s; i2++)
len += ab[i1][i2+r] * w1[i2];
len = sqrt(len);
for (i2 = 0; i2 < s; i2++)
ab[i1][i2+r] /= len;
}
}
}
else
sw = 0;
// 領域の解放
for (i1 = 0; i1 < q; i1++)
delete [] C[i1];
for (i1 = 0; i1 < r; i1++) {
delete [] B[i1];
delete [] C11[i1];
delete [] C11i[i1];
delete [] C12[i1];
}
for (i1 = 0; i1 < s; i1++) {
delete [] A[i1];
delete [] C21[i1];
delete [] C22[i1];
delete [] C22i[i1];
delete [] w[i1];
}
delete [] mean;
delete [] w1;
delete [] v1;
delete [] v2;
delete [] A;
delete [] B;
delete [] C;
delete [] C11;
delete [] C11i;
delete [] C12;
delete [] C21;
delete [] C22;
delete [] C22i;
delete [] w;
return sw;
}
/*******************************************************/
/* 線形連立方程式を解く（逆行列を求める） */
/* w : 方程式の左辺及び右辺 */
/* n : 方程式の数 */
/* m : 方程式の右辺の列の数 */
/* eps : 正則性を判定する規準 */
/* return : =0 : 正常 */
/* =1 : 逆行列が存在しない */
/*******************************************************/
int gauss(double **w, int n, int m, double eps)
{
double y1, y2;
int ind = 0, nm, m1, m2, i1, i2, i3;
nm = n + m;
for (i1 = 0; i1 < n && ind == 0; i1++) {
y1 = .0;
m1 = i1 + 1;
m2 = 0;
for (i2 = i1; i2 < n; i2++) {
y2 = fabs(w[i2][i1]);
if (y1 < y2) {
y1 = y2;
m2 = i2;
}
}
if (y1 < eps)
ind = 1;
else {
for (i2 = i1; i2 < nm; i2++) {
y1 = w[i1][i2];
w[i1][i2] = w[m2][i2];
w[m2][i2] = y1;
}
y1 = 1.0 / w[i1][i1];
for (i2 = m1; i2 < nm; i2++)
w[i1][i2] *= y1;
for (i2 = 0; i2 < n; i2++) {
if (i2 != i1) {
for (i3 = m1; i3 < nm; i3++)
w[i2][i3] -= w[i2][i1] * w[i1][i3];
}
}
}
}
return(ind);
}
/****************************************/
/* 最大固有値と固有ベクトル（べき乗法） */
/* i : 何番目の固有値かを示す */
/* n : 次数 */
/* m : 丸め誤差の修正回数 */
/* ct : 最大繰り返し回数 */
/* eps : 収束判定条件 */
/* A : 対象とする行列 */
/* B : nxnの行列，最初は，単位行列 */
/* C : 作業域，nxnの行列 */
/* r : 固有値 */
/* v : 各行が固有ベクトル（nxn） */
/* v0 : 固有ベクトルの初期値 */
/* v1,v2 : 作業域（n次元ベクトル） */
/* return : 求まった固有値の数 */
/* coded by Y.Suganuma */
/****************************************/
#include <math.h>
int power(int i, int n, int m, int ct, double eps, double **A, double **B,
double **C, double *r, double **v, double *v0, double *v1, double *v2)
{
double l1, l2, x;
int i1, i2, i3, ind, k, k1;
// 初期設定
ind = i;
k = 0;
l1 = 0.0;
for (i1 = 0; i1 < n; i1++) {
l1 += v0[i1] * v0[i1];
v1[i1] = v0[i1];
}
l1 = sqrt(l1);
// 繰り返し計算
while (k < ct) {
// 丸め誤差の修正
if (k%m == 0) {
l2 = 0.0;
for (i1 = 0; i1 < n; i1++) {
v2[i1] = 0.0;
for (i2 = 0; i2 < n; i2++)
v2[i1] += B[i1][i2] * v1[i2];
l2 += v2[i1] * v2[i1];
}
l2 = sqrt(l2);
for (i1 = 0; i1 < n; i1++)
v1[i1] = v2[i1] / l2;
}
// 次の近似
l2 = 0.0;
for (i1 = 0; i1 < n; i1++) {
v2[i1] = 0.0;
for (i2 = 0; i2 < n; i2++)
v2[i1] += A[i1][i2] * v1[i2];
l2 += v2[i1] * v2[i1];
}
l2 = sqrt(l2);
for (i1 = 0; i1 < n; i1++)
v2[i1] /= l2;
// 収束判定
// 収束した場合
if (fabs((l2-l1)/l1) < eps) {
k1 = -1;
for (i1 = 0; i1 < n && k1 < 0; i1++) {
if (fabs(v2[i1]) > 0.001) {
k1 = i1;
if (v2[k1]*v1[k1] < 0.0)
l2 = -l2;
}
}
k = ct;
r[i] = l2;
for (i1 = 0; i1 < n; i1++)
v[i][i1] = v2[i1];
if (i == n-1)
ind = i + 1;
else {
for (i1 = 0; i1 < n; i1++) {
for (i2 = 0; i2 < n; i2++) {
C[i1][i2] = 0.0;
for (i3 = 0; i3 < n; i3++) {
x = (i1 == i3) ? A[i1][i3] - l2 : A[i1][i3];
C[i1][i2] += x * B[i3][i2];
}
}
}
for (i1 = 0; i1 < n; i1++) {
for (i2 = 0; i2 < n; i2++)
B[i1][i2] = C[i1][i2];
}
for (i1 = 0; i1 < n; i1++) {
v1[i1] = 0.0;
for (i2 = 0; i2 < n; i2++)
v1[i1] += B[i1][i2] * v0[i2];
}
for (i1 = 0; i1 < n; i1++)
v0[i1] = v1[i1];
ind = power(i+1, n, m, ct, eps, A, B, C, r, v, v0, v1, v2);
}
}
// 収束しない場合
else {
for (i1 = 0; i1 < n; i1++)
v1[i1] = v2[i1];
l1 = l2;
k++;
}
}
return ind;
}
/*
---------データ例（コメント部分を除いて下さい）---------
2 2 100 // 各組の変数の数（r と s， r ≦ s）とデータの数（n）
66 22 44 31 // x1, x2, x3, x4
25 74 17 81
50 23 53 71
25 57 19 81
74 47 64 47
39 33 48 46
14 22 9 69
67 60 49 26
42 40 77 65
11 80 0 86
32 0 43 74
68 69 44 68
24 49 9 71
42 74 28 46
60 58 73 28
36 37 33 68
24 44 19 83
30 40 31 50
55 40 60 49
63 47 94 41
72 30 100 45
19 22 13 75
43 39 43 34
90 83 92 31
51 77 52 82
53 70 34 31
28 51 53 44
40 62 42 79
31 48 22 68
57 29 51 30
64 89 57 42
49 82 72 29
53 31 55 43
79 52 70 10
45 19 43 57
35 34 34 89
4 69 0 100
49 49 66 66
92 82 97 6
5 89 0 100
65 26 83 28
56 36 64 38
48 50 25 22
30 30 15 55
40 65 38 42
14 67 9 67
84 96 90 8
53 64 51 54
50 89 60 52
76 41 68 9
49 40 53 49
78 66 66 17
76 58 90 29
41 15 40 49
63 60 55 33
40 36 49 67
78 54 71 18
62 72 69 12
64 47 42 53
56 64 9 15
77 35 56 25
44 12 46 87
80 9 56 19
36 21 52 78
48 63 64 48
43 61 50 47
58 23 28 50
90 12 100 0
13 33 11 77
67 44 48 28
75 45 68 17
81 22 89 9
46 45 59 55
56 49 64 55
65 62 72 27
34 49 29 77
45 33 60 63
20 45 14 99
33 38 26 87
44 51 69 52
64 57 64 48
44 64 51 28
63 48 56 11
29 39 33 84
40 48 51 54
40 38 26 62
68 46 61 26
58 45 68 48
64 44 77 63
59 62 44 66
81 53 93 19
23 34 12 68
51 35 55 46
74 70 84 17
42 33 56 44
46 31 46 53
33 57 38 63
40 24 20 42
53 36 60 31
0 34 0 100
*/
/****************************/
/* 正準相関分析 */
/* coded by Y.Suganuma */
/****************************/
import java.io.*;
import java.util.StringTokenizer;
public class Test {
public static void main(String args[]) throws IOException
{
double x[][], ryz[], ab[][], v0[];
int i1, i2, q, r, s, n, sw;
StringTokenizer str;
BufferedReader in = new BufferedReader(new InputStreamReader(System.in));
// 各組の変数の数とデータの数
str = new StringTokenizer(in.readLine(), " ");
r = Integer.parseInt(str.nextToken());
s = Integer.parseInt(str.nextToken());
n = Integer.parseInt(str.nextToken());
q = r + s;
ryz = new double [r];
v0 = new double [r];
x = new double [q][n];
ab = new double [r][q];
for (i1 = 0; i1 < r; i1++)
v0[i1] = 1.0;
// データ
for (i1 = 0; i1 < n; i1++) {
str = new StringTokenizer(in.readLine(), " ");
for (i2 = 0; i2 < q; i2++)
x[i2][i1] = Double.parseDouble(str.nextToken());
}
sw = canonical(r, s, n, x, ryz, ab, 1.0e-10, v0, 15, 200);
if (sw > 0) {
for (i1 = 0; i1 < sw; i1++) {
System.out.println("相関係数 " + ryz[i1]);
System.out.print(" a");
for (i2 = 0; i2 < r; i2++)
System.out.print(" " + ab[i1][i2]);
System.out.println();
System.out.print(" b");
for (i2 = 0; i2 < s; i2++)
System.out.print(" " + ab[i1][r+i2]);
System.out.println();
}
}
else
System.out.println("***error 解を求めることができませんでした");
}
/***********************************/
/* 正準相関分析 */
/* r,s : 各組における変数の数 */
/* n : データの数 */
/* x : データ */
/* ryz : 相関係数 */
/* ab : 各組の係数（a,bの順） */
/* eps : 正則性を判定する規準 */
/* v0 : 固有ベクトルの初期値 */
/* m : 丸め誤差の修正回数 */
/* ct : 最大繰り返し回数 */
/* return : >0 : 相関係数の数 */
/* =0 : エラー */
/***********************************/
static int canonical(int r, int s, int n, double x[][], double ryz[],
double ab[][], double eps, double v0[], int m, int ct)
{
double A[][], B[][], C[][], C11[][], C11i[][], C12[][], C21[][], C22[][], C22i[][],
mean[], vv, w[][], w1[], len, v1[], v2[];
int i1, i2, i3, q, sw = 0, sw1;
// 領域の確保
q = r + s;
mean = new double [q];
w1 = new double [s];
v1 = new double [r];
v2 = new double [r];
A = new double [s][s];
B = new double [r][r];
C = new double [q][q];
C11 = new double [r][r];
C11i = new double [r][r];
C12 = new double [r][s];
C21 = new double [s][r];
C22 = new double [s][s];
C22i = new double [s][s];
w = new double [s][2*s];
// 平均値の計算
for (i1 = 0; i1 < q; i1++) {
mean[i1] = 0.0;
for (i2 = 0; i2 < n; i2++)
mean[i1] += x[i1][i2];
mean[i1] /= n;
}
// 分散強分散行列の計算
for (i1 = 0; i1 < q; i1++) {
for (i2 = i1; i2 < q; i2++) {
vv = 0.0;
for (i3 = 0; i3 < n; i3++)
vv += (x[i1][i3] - mean[i1]) * (x[i2][i3] - mean[i2]);
vv /= (n - 1);
C[i1][i2] = vv;
if (i1 != i2)
C[i2][i1] = vv;
}
}
// C11, C12, C21, C22 の設定
// C12
for (i1 = 0; i1 < r; i1++) {
for (i2 = 0; i2 < s; i2++)
C12[i1][i2] = C[i1][i2+r];
}
// C21
for (i1 = 0; i1 < s; i1++) {
for (i2 = 0; i2 < r; i2++)
C21[i1][i2] = C[i1+r][i2];
}
// C11とその逆行列
for (i1 = 0; i1 < r; i1++) {
for (i2 = 0; i2 < r; i2++) {
w[i1][i2] = C[i1][i2];
C11[i1][i2] = C[i1][i2];
}
for (i2 = r; i2 < 2*r; i2++)
w[i1][i2] = 0.0;
w[i1][i1+r] = 1.0;
}
sw1 = gauss(w, r, r, 1.0e-10);
if (sw1 == 0) {
for (i1 = 0; i1 < r; i1++) {
for (i2 = 0; i2 < r; i2++)
C11i[i1][i2] = w[i1][i2+r];
}
}
else
sw = 1;
// C22とその逆行列
for (i1 = 0; i1 < s; i1++) {
for (i2 = 0; i2 < s; i2++) {
w[i1][i2] = C[i1+r][i2+r];
C22[i1][i2] = C[i1+r][i2+r];
}
for (i2 = s; i2 < 2*s; i2++)
w[i1][i2] = 0.0;
w[i1][i1+s] = 1.0;
}
sw1 = gauss(w, s, s, eps);
if (sw1 == 0) {
for (i1 = 0; i1 < s; i1++) {
for (i2 = 0; i2 < s; i2++)
C22i[i1][i2] = w[i1][i2+s];
}
}
else
sw = 1;
// 固有値λ及び固有ベクトルaの計算
if (sw == 0) {
// 行列の計算
for (i1 = 0; i1 < s; i1++) {
for (i2 = 0; i2 < r; i2++) {
A[i1][i2] = 0.0;
for (i3 = 0; i3 < s; i3++)
A[i1][i2] += C22i[i1][i3] * C21[i3][i2];
}
}
for (i1 = 0; i1 < r; i1++) {
for (i2 = 0; i2 < s; i2++) {
w[i1][i2] = 0.0;
for (i3 = 0; i3 < s; i3++)
w[i1][i2] += C12[i1][i3] * A[i3][i2];
}
}
for (i1 = 0; i1 < r; i1++) {
for (i2 = 0; i2 < r; i2++) {
A[i1][i2] = 0.0;
for (i3 = 0; i3 < r; i3++)
A[i1][i2] += C11i[i1][i3] * w[i3][i2];
}
}
// 固有値と固有ベクトル（べき乗法）
for (i1 = 0; i1 < r; i1++) {
for (i2 = 0; i2 < r; i2++)
B[i1][i2] = 0.0;
B[i1][i1] = 1.0;
}
sw = power(0, r, m, ct, eps, A, B, C, ryz, ab, v0, v1, v2);
if (sw > 0) {
for (i1 = 0; i1 < r; i1++) {
// 相関係数
ryz[i1] = Math.sqrt(ryz[i1]);
// 大きさの調整(a)
for (i2 = 0; i2 < r; i2++) {
w1[i2] = 0.0;
for (i3 = 0; i3 < r; i3++)
w1[i2] += C11[i2][i3] * ab[i1][i3];
}
len = 0.0;
for (i2 = 0; i2 < r; i2++)
len += ab[i1][i2] * w1[i2];
len = Math.sqrt(len);
for (i2 = 0; i2 < r; i2++)
ab[i1][i2] /= len;
// bの計算
for (i2 = 0; i2 < s; i2++) {
w1[i2] = 0.0;
for (i3 = 0; i3 < r; i3++)
w1[i2] += C21[i2][i3] * ab[i1][i3];
}
for (i2 = 0; i2 < s; i2++) {
ab[i1][i2+r] = 0.0;
for (i3 = 0; i3 < s; i3++)
ab[i1][i2+r] += C22i[i2][i3] * w1[i3];
}
for (i2 = 0; i2 < s; i2++)
ab[i1][i2+r] /= ryz[i1];
// 大きさの調整(b)
for (i2 = 0; i2 < s; i2++) {
w1[i2] = 0.0;
for (i3 = 0; i3 < s; i3++)
w1[i2] += C22[i2][i3] * ab[i1][i3+r];
}
len = 0.0;
for (i2 = 0; i2 < s; i2++)
len += ab[i1][i2+r] * w1[i2];
len = Math.sqrt(len);
for (i2 = 0; i2 < s; i2++)
ab[i1][i2+r] /= len;
}
}
}
else
sw = 0;
return sw;
}
/****************************************/
/* 最大固有値と固有ベクトル（べき乗法） */
/* i : 何番目の固有値かを示す */
/* n : 次数 */
/* m : 丸め誤差の修正回数 */
/* ct : 最大繰り返し回数 */
/* eps : 収束判定条件 */
/* A : 対象とする行列 */
/* B : nxnの行列，最初は，単位行列 */
/* C : 作業域，nxnの行列 */
/* r : 固有値 */
/* v : 各行が固有ベクトル（nxn） */
/* v0 : 固有ベクトルの初期値 */
/* v1,v2 : 作業域（n次元ベクトル） */
/* return : 求まった固有値の数 */
/* coded by Y.Suganuma */
/****************************************/
static int power(int i, int n, int m, int ct, double eps, double A[][], double B[][],
double C[][], double r[], double v[][], double v0[], double v1[], double v2[])
{
double l1, l2, x;
int i1, i2, i3, ind, k, k1;
// 初期設定
ind = i;
k = 0;
l1 = 0.0;
for (i1 = 0; i1 < n; i1++) {
l1 += v0[i1] * v0[i1];
v1[i1] = v0[i1];
}
l1 = Math.sqrt(l1);
// 繰り返し計算
while (k < ct) {
// 丸め誤差の修正
if (k%m == 0) {
l2 = 0.0;
for (i1 = 0; i1 < n; i1++) {
v2[i1] = 0.0;
for (i2 = 0; i2 < n; i2++)
v2[i1] += B[i1][i2] * v1[i2];
l2 += v2[i1] * v2[i1];
}
l2 = Math.sqrt(l2);
for (i1 = 0; i1 < n; i1++)
v1[i1] = v2[i1] / l2;
}
// 次の近似
l2 = 0.0;
for (i1 = 0; i1 < n; i1++) {
v2[i1] = 0.0;
for (i2 = 0; i2 < n; i2++)
v2[i1] += A[i1][i2] * v1[i2];
l2 += v2[i1] * v2[i1];
}
l2 = Math.sqrt(l2);
for (i1 = 0; i1 < n; i1++)
v2[i1] /= l2;
// 収束判定
// 収束した場合
if (Math.abs((l2-l1)/l1) < eps) {
k1 = -1;
for (i1 = 0; i1 < n && k1 < 0; i1++) {
if (Math.abs(v2[i1]) > 0.001) {
k1 = i1;
if (v2[k1]*v1[k1] < 0.0)
l2 = -l2;
}
}
k = ct;
r[i] = l2;
for (i1 = 0; i1 < n; i1++)
v[i][i1] = v2[i1];
if (i == n-1)
ind = i + 1;
else {
for (i1 = 0; i1 < n; i1++) {
for (i2 = 0; i2 < n; i2++) {
C[i1][i2] = 0.0;
for (i3 = 0; i3 < n; i3++) {
x = (i1 == i3) ? A[i1][i3] - l2 : A[i1][i3];
C[i1][i2] += x * B[i3][i2];
}
}
}
for (i1 = 0; i1 < n; i1++) {
for (i2 = 0; i2 < n; i2++)
B[i1][i2] = C[i1][i2];
}
for (i1 = 0; i1 < n; i1++) {
v1[i1] = 0.0;
for (i2 = 0; i2 < n; i2++)
v1[i1] += B[i1][i2] * v0[i2];
}
for (i1 = 0; i1 < n; i1++)
v0[i1] = v1[i1];
ind = power(i+1, n, m, ct, eps, A, B, C, r, v, v0, v1, v2);
}
}
// 収束しない場合
else {
for (i1 = 0; i1 < n; i1++)
v1[i1] = v2[i1];
l1 = l2;
k++;
}
}
return ind;
}
/*******************************************************/
/* 線形連立方程式を解く（逆行列を求める） */
/* w : 方程式の左辺及び右辺 */
/* n : 方程式の数 */
/* m : 方程式の右辺の列の数 */
/* eps : 逆行列の存在を判定する規準 */
/* return : =0 : 正常 */
/* =1 : 逆行列が存在しない */
/*******************************************************/
static int gauss(double w[][], int n, int m, double eps)
{
double y1, y2;
int ind = 0, nm, m1, m2, i1, i2, i3;
nm = n + m;
for (i1 = 0; i1 < n && ind == 0; i1++) {
y1 = .0;
m1 = i1 + 1;
m2 = 0;
// ピボット要素の選択
for (i2 = i1; i2 < n; i2++) {
y2 = Math.abs(w[i2][i1]);
if (y1 < y2) {
y1 = y2;
m2 = i2;
}
}
// 逆行列が存在しない
if (y1 < eps)
ind = 1;
// 逆行列が存在する
else {
// 行の入れ替え
for (i2 = i1; i2 < nm; i2++) {
y1 = w[i1][i2];
w[i1][i2] = w[m2][i2];
w[m2][i2] = y1;
}
// 掃き出し操作
y1 = 1.0 / w[i1][i1];
for (i2 = m1; i2 < nm; i2++)
w[i1][i2] *= y1;
for (i2 = 0; i2 < n; i2++) {
if (i2 != i1) {
for (i3 = m1; i3 < nm; i3++)
w[i2][i3] -= w[i2][i1] * w[i1][i3];
}
}
}
}
return(ind);
}
}
/*
---------データ例（コメント部分を除いて下さい）---------
2 2 100 // 各組の変数の数（r と s， r ≦ s）とデータの数（n）
66 22 44 31 // x1, x2, x3, x4
25 74 17 81
50 23 53 71
25 57 19 81
74 47 64 47
39 33 48 46
14 22 9 69
67 60 49 26
42 40 77 65
11 80 0 86
32 0 43 74
68 69 44 68
24 49 9 71
42 74 28 46
60 58 73 28
36 37 33 68
24 44 19 83
30 40 31 50
55 40 60 49
63 47 94 41
72 30 100 45
19 22 13 75
43 39 43 34
90 83 92 31
51 77 52 82
53 70 34 31
28 51 53 44
40 62 42 79
31 48 22 68
57 29 51 30
64 89 57 42
49 82 72 29
53 31 55 43
79 52 70 10
45 19 43 57
35 34 34 89
4 69 0 100
49 49 66 66
92 82 97 6
5 89 0 100
65 26 83 28
56 36 64 38
48 50 25 22
30 30 15 55
40 65 38 42
14 67 9 67
84 96 90 8
53 64 51 54
50 89 60 52
76 41 68 9
49 40 53 49
78 66 66 17
76 58 90 29
41 15 40 49
63 60 55 33
40 36 49 67
78 54 71 18
62 72 69 12
64 47 42 53
56 64 9 15
77 35 56 25
44 12 46 87
80 9 56 19
36 21 52 78
48 63 64 48
43 61 50 47
58 23 28 50
90 12 100 0
13 33 11 77
67 44 48 28
75 45 68 17
81 22 89 9
46 45 59 55
56 49 64 55
65 62 72 27
34 49 29 77
45 33 60 63
20 45 14 99
33 38 26 87
44 51 69 52
64 57 64 48
44 64 51 28
63 48 56 11
29 39 33 84
40 48 51 54
40 38 26 62
68 46 61 26
58 45 68 48
64 44 77 63
59 62 44 66
81 53 93 19
23 34 12 68
51 35 55 46
74 70 84 17
42 33 56 44
46 31 46 53
33 57 38 63
40 24 20 42
53 36 60 31
0 34 0 100
*/
<!DOCTYPE HTML>
<HTML>
<HEAD>
<TITLE>正準相関分析</TITLE>
<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=utf-8">
<SCRIPT TYPE="text/javascript">
function main()
{
// データ
let r = parseInt(document.getElementById("r").value);
let s = parseInt(document.getElementById("s").value);
let n = parseInt(document.getElementById("n").value);
let q = r + s;
let ryz = new Array();
let v0 = new Array();
let ab = new Array();
for (let i1 = 0; i1 < r; i1++) {
v0[i1] = 1.0;
ab[i1] = new Array();
}
let x = new Array();
for (let i1 = 0; i1 < q; i1++)
x[i1] = new Array();
let ss = (document.getElementById("data").value).split(/ {1,}|\n{1,}/);
let k = 0;
for (let i1 = 0; i1 < n; i1++) {
for (let i2 = 0; i2 < q; i2++) {
x[i2][i1] = parseFloat(ss[k]);
k++;
}
}
// 計算
let sw = canonical(r, s, n, x, ryz, ab, 1.0e-10, v0, 15, 200);
if (sw == 0)
document.getElementById("ans").value = "逆行列が存在しない";
else {
let str = "";
for (let i1 = 0; i1 < sw; i1++) {
str = str + "相関係数 " + ryz[i1] + "\n";
str = str + " a";
for (let i2 = 0; i2 < r; i2++)
str = str + " " + ab[i1][i2];
str = str + "\n b";
for (let i2 = 0; i2 < s; i2++)
str = str + " " + ab[i1][r+i2];
str = str + "\n";
}
document.getElementById("ans").value = str;
}
}
/***********************************/
/* 正準相関分析 */
/* r,s : 各組における変数の数 */
/* n : データの数 */
/* x : データ */
/* ryz : 相関係数 */
/* ab : 各組の係数（a,bの順） */
/* eps : 正則性を判定する規準 */
/* v0 : 固有ベクトルの初期値 */
/* m : 丸め誤差の修正回数 */
/* ct : 最大繰り返し回数 */
/* return : >0 : 相関係数の数 */
/* =0 : エラー */
/***********************************/
function canonical(r, s, n, x, ryz, ab, eps, v0, m, ct)
{
let vv;
let len;
let i1;
let i2;
let i3;
let q = r + s;
let sw = 0;
let sw1;
// 領域の確保
let mean = new Array();
let w1 = new Array();
let v1 = new Array();
let v2 = new Array();
let C = new Array();
for (i1 = 0; i1 < q; i1++)
C[i1] = new Array();
let B = new Array();
let C11 = new Array();
let C11i = new Array();
let C12 = new Array();
for (i1 = 0; i1 < r; i1++) {
B[i1] = new Array();
C11[i1] = new Array();
C11i[i1] = new Array();
C12[i1] = new Array();
}
let A = new Array();
let C21 = new Array();
let C22 = new Array();
let C22i = new Array();
let w = new Array();
for (i1 = 0; i1 < s; i1++) {
A[i1] = new Array();
C21[i1] = new Array();
C22[i1] = new Array();
C22i[i1] = new Array();
w[i1] = new Array();
}
// 平均値の計算
for (i1 = 0; i1 < q; i1++) {
mean[i1] = 0.0;
for (i2 = 0; i2 < n; i2++)
mean[i1] += x[i1][i2];
mean[i1] /= n;
}
// 分散強分散行列の計算
for (i1 = 0; i1 < q; i1++) {
for (i2 = i1; i2 < q; i2++) {
vv = 0.0;
for (i3 = 0; i3 < n; i3++)
vv += (x[i1][i3] - mean[i1]) * (x[i2][i3] - mean[i2]);
vv /= (n - 1);
C[i1][i2] = vv;
if (i1 != i2)
C[i2][i1] = vv;
}
}
// C11, C12, C21, C22 の設定
// C12
for (i1 = 0; i1 < r; i1++) {
for (i2 = 0; i2 < s; i2++)
C12[i1][i2] = C[i1][i2+r];
}
// C21
for (i1 = 0; i1 < s; i1++) {
for (i2 = 0; i2 < r; i2++)
C21[i1][i2] = C[i1+r][i2];
}
// C11とその逆行列
for (i1 = 0; i1 < r; i1++) {
for (i2 = 0; i2 < r; i2++) {
w[i1][i2] = C[i1][i2];
C11[i1][i2] = C[i1][i2];
}
for (i2 = r; i2 < 2*r; i2++)
w[i1][i2] = 0.0;
w[i1][i1+r] = 1.0;
}
sw1 = gauss(w, r, r, 1.0e-10);
if (sw1 == 0) {
for (i1 = 0; i1 < r; i1++) {
for (i2 = 0; i2 < r; i2++)
C11i[i1][i2] = w[i1][i2+r];
}
}
else
sw = 1;
// C22とその逆行列
for (i1 = 0; i1 < s; i1++) {
for (i2 = 0; i2 < s; i2++) {
w[i1][i2] = C[i1+r][i2+r];
C22[i1][i2] = C[i1+r][i2+r];
}
for (i2 = s; i2 < 2*s; i2++)
w[i1][i2] = 0.0;
w[i1][i1+s] = 1.0;
}
sw1 = gauss(w, s, s, eps);
if (sw1 == 0) {
for (i1 = 0; i1 < s; i1++) {
for (i2 = 0; i2 < s; i2++)
C22i[i1][i2] = w[i1][i2+s];
}
}
else
sw = 1;
// 固有値λ及び固有ベクトルaの計算
if (sw == 0) {
// 行列の計算
for (i1 = 0; i1 < s; i1++) {
for (i2 = 0; i2 < r; i2++) {
A[i1][i2] = 0.0;
for (i3 = 0; i3 < s; i3++)
A[i1][i2] += C22i[i1][i3] * C21[i3][i2];
}
}
for (i1 = 0; i1 < r; i1++) {
for (i2 = 0; i2 < s; i2++) {
w[i1][i2] = 0.0;
for (i3 = 0; i3 < s; i3++)
w[i1][i2] += C12[i1][i3] * A[i3][i2];
}
}
for (i1 = 0; i1 < r; i1++) {
for (i2 = 0; i2 < r; i2++) {
A[i1][i2] = 0.0;
for (i3 = 0; i3 < r; i3++)
A[i1][i2] += C11i[i1][i3] * w[i3][i2];
}
}
// 固有値と固有ベクトル（べき乗法）
for (i1 = 0; i1 < r; i1++) {
for (i2 = 0; i2 < r; i2++)
B[i1][i2] = 0.0;
B[i1][i1] = 1.0;
}
sw = power(0, r, m, ct, eps, A, B, C, ryz, ab, v0, v1, v2);
if (sw > 0) {
for (i1 = 0; i1 < r; i1++) {
// 相関係数
ryz[i1] = Math.sqrt(ryz[i1]);
// 大きさの調整(a)
for (i2 = 0; i2 < r; i2++) {
w1[i2] = 0.0;
for (i3 = 0; i3 < r; i3++)
w1[i2] += C11[i2][i3] * ab[i1][i3];
}
len = 0.0;
for (i2 = 0; i2 < r; i2++)
len += ab[i1][i2] * w1[i2];
len = Math.sqrt(len);
for (i2 = 0; i2 < r; i2++)
ab[i1][i2] /= len;
// bの計算
for (i2 = 0; i2 < s; i2++) {
w1[i2] = 0.0;
for (i3 = 0; i3 < r; i3++)
w1[i2] += C21[i2][i3] * ab[i1][i3];
}
for (i2 = 0; i2 < s; i2++) {
ab[i1][i2+r] = 0.0;
for (i3 = 0; i3 < s; i3++)
ab[i1][i2+r] += C22i[i2][i3] * w1[i3];
}
for (i2 = 0; i2 < s; i2++)
ab[i1][i2+r] /= ryz[i1];
// 大きさの調整(b)
for (i2 = 0; i2 < s; i2++) {
w1[i2] = 0.0;
for (i3 = 0; i3 < s; i3++)
w1[i2] += C22[i2][i3] * ab[i1][i3+r];
}
len = 0.0;
for (i2 = 0; i2 < s; i2++)
len += ab[i1][i2+r] * w1[i2];
len = Math.sqrt(len);
for (i2 = 0; i2 < s; i2++)
ab[i1][i2+r] /= len;
}
}
}
else
sw = 0;
return sw;
}
/****************************************/
/* 最大固有値と固有ベクトル（べき乗法） */
/* i : 何番目の固有値かを示す */
/* n : 次数 */
/* m : 丸め誤差の修正回数 */
/* ct : 最大繰り返し回数 */
/* eps : 収束判定条件 */
/* A : 対象とする行列 */
/* B : nxnの行列，最初は，単位行列 */
/* C : 作業域，nxnの行列 */
/* r : 固有値 */
/* v : 各行が固有ベクトル（nxn） */
/* v0 : 固有ベクトルの初期値 */
/* v1,v2 : 作業域（n次元ベクトル） */
/* return : 求まった固有値の数 */
/* coded by Y.Suganuma */
/****************************************/
function power(i, n, m, ct, eps, A, B, C, r, v, v0, v1, v2)
{
let l1;
let l2;
let x;
let i1;
let i2;
let i3;
let ind;
let k;
let k1;
// 初期設定
ind = i;
k = 0;
l1 = 0.0;
for (i1 = 0; i1 < n; i1++) {
l1 += v0[i1] * v0[i1];
v1[i1] = v0[i1];
}
l1 = Math.sqrt(l1);
// 繰り返し計算
while (k < ct) {
// 丸め誤差の修正
if (k%m == 0) {
l2 = 0.0;
for (i1 = 0; i1 < n; i1++) {
v2[i1] = 0.0;
for (i2 = 0; i2 < n; i2++)
v2[i1] += B[i1][i2] * v1[i2];
l2 += v2[i1] * v2[i1];
}
l2 = Math.sqrt(l2);
for (i1 = 0; i1 < n; i1++)
v1[i1] = v2[i1] / l2;
}
// 次の近似
l2 = 0.0;
for (i1 = 0; i1 < n; i1++) {
v2[i1] = 0.0;
for (i2 = 0; i2 < n; i2++)
v2[i1] += A[i1][i2] * v1[i2];
l2 += v2[i1] * v2[i1];
}
l2 = Math.sqrt(l2);
for (i1 = 0; i1 < n; i1++)
v2[i1] /= l2;
// 収束判定
// 収束した場合
if (Math.abs((l2-l1)/l1) < eps) {
k1 = -1;
for (i1 = 0; i1 < n && k1 < 0; i1++) {
if (Math.abs(v2[i1]) > 0.001) {
k1 = i1;
if (v2[k1]*v1[k1] < 0.0)
l2 = -l2;
}
}
k = ct;
r[i] = l2;
for (i1 = 0; i1 < n; i1++)
v[i][i1] = v2[i1];
if (i == n-1)
ind = i + 1;
else {
for (i1 = 0; i1 < n; i1++) {
for (i2 = 0; i2 < n; i2++) {
C[i1][i2] = 0.0;
for (i3 = 0; i3 < n; i3++) {
x = (i1 == i3) ? A[i1][i3] - l2 : A[i1][i3];
C[i1][i2] += x * B[i3][i2];
}
}
}
for (i1 = 0; i1 < n; i1++) {
for (i2 = 0; i2 < n; i2++)
B[i1][i2] = C[i1][i2];
}
for (i1 = 0; i1 < n; i1++) {
v1[i1] = 0.0;
for (i2 = 0; i2 < n; i2++)
v1[i1] += B[i1][i2] * v0[i2];
}
for (i1 = 0; i1 < n; i1++)
v0[i1] = v1[i1];
ind = power(i+1, n, m, ct, eps, A, B, C, r, v, v0, v1, v2);
}
}
// 収束しない場合
else {
for (i1 = 0; i1 < n; i1++)
v1[i1] = v2[i1];
l1 = l2;
k++;
}
}
return ind;
}
/******************************************/
/* 線形連立方程式を解く（逆行列を求める） */
/* w : 方程式の左辺及び右辺 */
/* n : 方程式の数 */
/* m : 方程式の右辺の列の数 */
/* eps : 逆行列の存在を判定する規準 */
/* return : =0 : 正常 */
/* =1 : 逆行列が存在しない */
/******************************************/
function gauss(w, n, m, eps)
{
let y1, y2;
let ind = 0;
let nm;
let m1;
let m2;
let i1;
let i2;
let i3;
nm = n + m;
for (i1 = 0; i1 < n && ind == 0; i1++) {
y1 = .0;
m1 = i1 + 1;
m2 = 0;
// ピボット要素の選択
for (i2 = i1; i2 < n; i2++) {
y2 = Math.abs(w[i2][i1]);
if (y1 < y2) {
y1 = y2;
m2 = i2;
}
}
// 逆行列が存在しない
if (y1 < eps)
ind = 1;
// 逆行列が存在する
else {
// 行の入れ替え
for (i2 = i1; i2 < nm; i2++) {
y1 = w[i1][i2];
w[i1][i2] = w[m2][i2];
w[m2][i2] = y1;
}
// 掃き出し操作
y1 = 1.0 / w[i1][i1];
for (i2 = m1; i2 < nm; i2++)
w[i1][i2] *= y1;
for (i2 = 0; i2 < n; i2++) {
if (i2 != i1) {
for (i3 = m1; i3 < nm; i3++)
w[i2][i3] -= w[i2][i1] * w[i1][i3];
}
}
}
}
return(ind);
}
</SCRIPT>
</HEAD>
<BODY STYLE="font-size: 130%; background-color: #eeffee;">
<H2 STYLE="text-align:center"><B>正準相関分析</B></H2>
<DL>
<DT> テキストフィールドおよびテキストエリアには，例として，テキストエリアに与えられた 100 個のデータに対して正準相関分析を行う場合に対する値が設定されています．他の問題を実行する場合は，それらを適切に修正してください．
</DL>
<DIV STYLE="text-align:center">
変数の数(r)：<INPUT ID="r" STYLE="font-size: 100%" TYPE="text" SIZE="2" VALUE="2">
変数の数(s)：<INPUT ID="s" STYLE="font-size: 100%" TYPE="text" SIZE="2" VALUE="2">
データの数：<INPUT ID="n" STYLE="font-size: 100%" TYPE="text" SIZE="2" VALUE="100"><BR><BR>
データ：<TEXTAREA ID="data" COLS="30" ROWS="15" STYLE="font-size: 100%">
66 22 44 31
25 74 17 81
50 23 53 71
25 57 19 81
74 47 64 47
39 33 48 46
14 22 9 69
67 60 49 26
42 40 77 65
11 80 0 86
32 0 43 74
68 69 44 68
24 49 9 71
42 74 28 46
60 58 73 28
36 37 33 68
24 44 19 83
30 40 31 50
55 40 60 49
63 47 94 41
72 30 100 45
19 22 13 75
43 39 43 34
90 83 92 31
51 77 52 82
53 70 34 31
28 51 53 44
40 62 42 79
31 48 22 68
57 29 51 30
64 89 57 42
49 82 72 29
53 31 55 43
79 52 70 10
45 19 43 57
35 34 34 89
4 69 0 100
49 49 66 66
92 82 97 6
5 89 0 100
65 26 83 28
56 36 64 38
48 50 25 22
30 30 15 55
40 65 38 42
14 67 9 67
84 96 90 8
53 64 51 54
50 89 60 52
76 41 68 9
49 40 53 49
78 66 66 17
76 58 90 29
41 15 40 49
63 60 55 33
40 36 49 67
78 54 71 18
62 72 69 12
64 47 42 53
56 64 9 15
77 35 56 25
44 12 46 87
80 9 56 19
36 21 52 78
48 63 64 48
43 61 50 47
58 23 28 50
90 12 100 0
13 33 11 77
67 44 48 28
75 45 68 17
81 22 89 9
46 45 59 55
56 49 64 55
65 62 72 27
34 49 29 77
45 33 60 63
20 45 14 99
33 38 26 87
44 51 69 52
64 57 64 48
44 64 51 28
63 48 56 11
29 39 33 84
40 48 51 54
40 38 26 62
68 46 61 26
58 45 68 48
64 44 77 63
59 62 44 66
81 53 93 19
23 34 12 68
51 35 55 46
74 70 84 17
42 33 56 44
46 31 46 53
33 57 38 63
40 24 20 42
53 36 60 31
0 34 0 100
</TEXTAREA>
<BUTTON STYLE="font-size: 100%; background-color: pink" onClick="main()">OK</BUTTON><BR><BR>
<TEXTAREA ID="ans" COLS="50" ROWS="5" STYLE="font-size: 100%;"></TEXTAREA>
</DIV>
</BODY>
</HTML>
<?php
/****************************/
/* 正準相関分析 */
/* coded by Y.Suganuma */
/****************************/
fscanf(STDIN, "%d %d %d", $r, $s, $n); // 各組の変数の数とデータの数
$q = $r + $s;
$ryz = array($r);
$v0 = array($r);
$x = array($q);
$ab = array($r);
for ($i1 = 0; $i1 < $q; $i1++)
$x[$i1] = array($n);
for ($i1 = 0; $i1 < $r; $i1++) {
$ab[$i1] = array($q);
$v0[$i1] = 1.0;
}
for ($i1 = 0; $i1 < $n; $i1++) { // データ
$str = trim(fgets(STDIN));
$x[0][$i1] = floatval(strtok($str, " "));
for ($i2 = 1; $i2 < $q; $i2++)
$x[$i2][$i1] = floatval(strtok(" "));
}
$sw = canonical($r, $s, $n, $x, $ryz, $ab, 1.0e-10, $v0, 15, 200);
if ($sw > 0) {
for ($i1 = 0; $i1 < $sw; $i1++) {
printf("相関係数 %f\n", $ryz[$i1]);
printf(" a");
for ($i2 = 0; $i2 < $r; $i2++)
printf(" %f", $ab[$i1][$i2]);
printf("\n");
printf(" b");
for ($i2 = 0; $i2 < $s; $i2++)
printf(" %f", $ab[$i1][$r+$i2]);
printf("\n");
}
}
else
printf("***error 解を求めることができませんでした\n");
/***********************************/
/* 正準相関分析 */
/* r,s : 各組における変数の数 */
/* n : データの数 */
/* x : データ */
/* ryz : 相関係数 */
/* ab : 各組の係数（a,bの順） */
/* eps : 正則性を判定する規準 */
/* v0 : 固有ベクトルの初期値 */
/* m : 丸め誤差の修正回数 */
/* ct : 最大繰り返し回数 */
/* return : >0 : 相関係数の数 */
/* =0 : エラー */
/***********************************/
function canonical($r, $s, $n, $x, &$ryz, &$ab, $eps, $v0, $m, $ct)
{
$sw = 0;
// 領域の確保
$q = $r + $s;
$mean = array($q);
$w1 = array($s);
$v1 = array($r);
$v2 = array($r);
$A = array($s);
$B = array($r);
$C = array($q);
$C11 = array($r);
$C11i = array($r);
$C12 = array($r);
$C21 = array($s);
$C22 = array($s);
$C22i = array($s);
$w = array($s);
for ($i1 = 0; $i1 < $q; $i1++)
$C[$i1] = array($q);
for ($i1 = 0; $i1 < $r; $i1++) {
$B[$i1] = array($r);
$C11[$i1] = array($r);
$C11i[$i1] = array($r);
$C12[$i1] = array($s);
}
for ($i1 = 0; $i1 < $s; $i1++) {
$A[$i1] = array($s);
$C21[$i1] = array($r);
$C22[$i1] = array($s);
$C22i[$i1] = array($s);
$w[$i1] = array(2*$s);
}
// 平均値の計算
for ($i1 = 0; $i1 < $q; $i1++) {
$mean[$i1] = 0.0;
for ($i2 = 0; $i2 < $n; $i2++)
$mean[$i1] += $x[$i1][$i2];
$mean[$i1] /= $n;
}
// 分散強分散行列の計算
for ($i1 = 0; $i1 < $q; $i1++) {
for ($i2 = $i1; $i2 < $q; $i2++) {
$vv = 0.0;
for ($i3 = 0; $i3 < $n; $i3++)
$vv += ($x[$i1][$i3] - $mean[$i1]) * ($x[$i2][$i3] - $mean[$i2]);
$vv /= ($n - 1);
$C[$i1][$i2] = $vv;
if ($i1 != $i2)
$C[$i2][$i1] = $vv;
}
}
// C11, C12, C21, C22 の設定
// C12
for ($i1 = 0; $i1 < $r; $i1++) {
for ($i2 = 0; $i2 < $s; $i2++)
$C12[$i1][$i2] = $C[$i1][$i2+$r];
}
// C21
for ($i1 = 0; $i1 < $s; $i1++) {
for ($i2 = 0; $i2 < $r; $i2++)
$C21[$i1][$i2] = $C[$i1+$r][$i2];
}
// C11とその逆行列
for ($i1 = 0; $i1 < $r; $i1++) {
for ($i2 = 0; $i2 < $r; $i2++) {
$w[$i1][$i2] = $C[$i1][$i2];
$C11[$i1][$i2] = $C[$i1][$i2];
}
for ($i2 = $r; $i2 < 2*$r; $i2++)
$w[$i1][$i2] = 0.0;
$w[$i1][$i1+$r] = 1.0;
}
$sw1 = gauss($w, $r, $r, 1.0e-10);
if ($sw1 == 0) {
for ($i1 = 0; $i1 < $r; $i1++) {
for ($i2 = 0; $i2 < $r; $i2++)
$C11i[$i1][$i2] = $w[$i1][$i2+$r];
}
}
else
$sw = 1;
// C22とその逆行列
for ($i1 = 0; $i1 < $s; $i1++) {
for ($i2 = 0; $i2 < $s; $i2++) {
$w[$i1][$i2] = $C[$i1+$r][$i2+$r];
$C22[$i1][$i2] = $C[$i1+$r][$i2+$r];
}
for ($i2 = $s; $i2 < 2*$s; $i2++)
$w[$i1][$i2] = 0.0;
$w[$i1][$i1+$s] = 1.0;
}
$sw1 = gauss($w, $s, $s, $eps);
if ($sw1 == 0) {
for ($i1 = 0; $i1 < $s; $i1++) {
for ($i2 = 0; $i2 < $s; $i2++)
$C22i[$i1][$i2] = $w[$i1][$i2+$s];
}
}
else
$sw = 1;
// 固有値λ及び固有ベクトルaの計算
if ($sw == 0) {
// 行列の計算
for ($i1 = 0; $i1 < $s; $i1++) {
for ($i2 = 0; $i2 < $r; $i2++) {
$A[$i1][$i2] = 0.0;
for ($i3 = 0; $i3 < $s; $i3++)
$A[$i1][$i2] += $C22i[$i1][$i3] * $C21[$i3][$i2];
}
}
for ($i1 = 0; $i1 < $r; $i1++) {
for ($i2 = 0; $i2 < $s; $i2++) {
$w[$i1][$i2] = 0.0;
for ($i3 = 0; $i3 < $s; $i3++)
$w[$i1][$i2] += $C12[$i1][$i3] * $A[$i3][$i2];
}
}
for ($i1 = 0; $i1 < $r; $i1++) {
for ($i2 = 0; $i2 < $r; $i2++) {
$A[$i1][$i2] = 0.0;
for ($i3 = 0; $i3 < $r; $i3++)
$A[$i1][$i2] += $C11i[$i1][$i3] * $w[$i3][$i2];
}
}
// 固有値と固有ベクトル（べき乗法）
for ($i1 = 0; $i1 < $r; $i1++) {
for ($i2 = 0; $i2 < $r; $i2++)
$B[$i1][$i2] = 0.0;
$B[$i1][$i1] = 1.0;
}
$sw = power(0, $r, $m, $ct, $eps, $A, $B, $C, $ryz, $ab, $v0, $v1, $v2);
if ($sw > 0) {
for ($i1 = 0; $i1 < $r; $i1++) {
// 相関係数
$ryz[$i1] = sqrt($ryz[$i1]);
// 大きさの調整(a)
for ($i2 = 0; $i2 < $r; $i2++) {
$w1[$i2] = 0.0;
for ($i3 = 0; $i3 < $r; $i3++)
$w1[$i2] += $C11[$i2][$i3] * $ab[$i1][$i3];
}
$len = 0.0;
for ($i2 = 0; $i2 < $r; $i2++)
$len += $ab[$i1][$i2] * $w1[$i2];
$len = sqrt($len);
for ($i2 = 0; $i2 < $r; $i2++)
$ab[$i1][$i2] /= $len;
// bの計算
for ($i2 = 0; $i2 < $s; $i2++) {
$w1[$i2] = 0.0;
for ($i3 = 0; $i3 < $r; $i3++)
$w1[$i2] += $C21[$i2][$i3] * $ab[$i1][$i3];
}
for ($i2 = 0; $i2 < $s; $i2++) {
$ab[$i1][$i2+$r] = 0.0;
for ($i3 = 0; $i3 < $s; $i3++)
$ab[$i1][$i2+$r] += $C22i[$i2][$i3] * $w1[$i3];
}
for ($i2 = 0; $i2 < $s; $i2++)
$ab[$i1][$i2+$r] /= $ryz[$i1];
// 大きさの調整(b)
for ($i2 = 0; $i2 < $s; $i2++) {
$w1[$i2] = 0.0;
for ($i3 = 0; $i3 < $s; $i3++)
$w1[$i2] += $C22[$i2][$i3] * $ab[$i1][$i3+$r];
}
$len = 0.0;
for ($i2 = 0; $i2 < $s; $i2++)
$len += $ab[$i1][$i2+$r] * $w1[$i2];
$len = sqrt($len);
for ($i2 = 0; $i2 < $s; $i2++)
$ab[$i1][$i2+$r] /= $len;
}
}
}
else
$sw = 0;
return $sw;
}
/*******************************************************/
/* 線形連立方程式を解く（逆行列を求める） */
/* w : 方程式の左辺及び右辺 */
/* n : 方程式の数 */
/* m : 方程式の右辺の列の数 */
/* eps : 正則性を判定する規準 */
/* return : =0 : 正常 */
/* =1 : 逆行列が存在しない */
/*******************************************************/
function gauss(&$w, $n, $m, $eps)
{
$ind = 0;
$nm = $n + $m;
for ($i1 = 0; $i1 < $n && $ind == 0; $i1++) {
$y1 = .0;
$m1 = $i1 + 1;
$m2 = 0;
for ($i2 = $i1; $i2 < $n; $i2++) {
$y2 = abs($w[$i2][$i1]);
if ($y1 < $y2) {
$y1 = $y2;
$m2 = $i2;
}
}
if ($y1 < $eps)
$ind = 1;
else {
for ($i2 = $i1; $i2 < $nm; $i2++) {
$y1 = $w[$i1][$i2];
$w[$i1][$i2] = $w[$m2][$i2];
$w[$m2][$i2] = $y1;
}
$y1 = 1.0 / $w[$i1][$i1];
for ($i2 = $m1; $i2 < $nm; $i2++)
$w[$i1][$i2] *= $y1;
for ($i2 = 0; $i2 < $n; $i2++) {
if ($i2 != $i1) {
for ($i3 = $m1; $i3 < $nm; $i3++)
$w[$i2][$i3] -= $w[$i2][$i1] * $w[$i1][$i3];
}
}
}
}
return($ind);
}
/****************************************/
/* 最大固有値と固有ベクトル（べき乗法） */
/* i : 何番目の固有値かを示す */
/* n : 次数 */
/* m : 丸め誤差の修正回数 */
/* ct : 最大繰り返し回数 */
/* eps : 収束判定条件 */
/* A : 対象とする行列 */
/* B : nxnの行列，最初は，単位行列 */
/* C : 作業域，nxnの行列 */
/* r : 固有値 */
/* v : 各行が固有ベクトル（nxn） */
/* v0 : 固有ベクトルの初期値 */
/* v1,v2 : 作業域（n次元ベクトル） */
/* return : 求まった固有値の数 */
/* coded by Y.Suganuma */
/****************************************/
function power($i, $n, $m, $ct, $eps, $A, $B, $C, &$r, &$v, $v0, $v1, $v2)
{
// 初期設定
$ind = $i;
$k = 0;
$l1 = 0.0;
for ($i1 = 0; $i1 < $n; $i1++) {
$l1 += $v0[$i1] * $v0[$i1];
$v1[$i1] = $v0[$i1];
}
$l1 = sqrt($l1);
// 繰り返し計算
while ($k < $ct) {
// 丸め誤差の修正
if ($k%$m == 0) {
$l2 = 0.0;
for ($i1 = 0; $i1 < $n; $i1++) {
$v2[$i1] = 0.0;
for ($i2 = 0; $i2 < $n; $i2++)
$v2[$i1] += $B[$i1][$i2] * $v1[$i2];
$l2 += $v2[$i1] * $v2[$i1];
}
$l2 = sqrt($l2);
for ($i1 = 0; $i1 < $n; $i1++)
$v1[$i1] = $v2[$i1] / $l2;
}
// 次の近似
$l2 = 0.0;
for ($i1 = 0; $i1 < $n; $i1++) {
$v2[$i1] = 0.0;
for ($i2 = 0; $i2 < $n; $i2++)
$v2[$i1] += $A[$i1][$i2] * $v1[$i2];
$l2 += $v2[$i1] * $v2[$i1];
}
$l2 = sqrt($l2);
for ($i1 = 0; $i1 < $n; $i1++)
$v2[$i1] /= $l2;
// 収束判定
// 収束した場合
if (abs(($l2-$l1)/$l1) < $eps) {
$k1 = -1;
for ($i1 = 0; $i1 < $n && $k1 < 0; $i1++) {
if (abs($v2[$i1]) > 0.001) {
$k1 = $i1;
if ($v2[$k1]*$v1[$k1] < 0.0)
$l2 = -$l2;
}
}
$k = $ct;
$r[$i] = $l2;
for ($i1 = 0; $i1 < $n; $i1++)
$v[$i][$i1] = $v2[$i1];
if ($i == $n-1)
$ind = $i + 1;
else {
for ($i1 = 0; $i1 < $n; $i1++) {
for ($i2 = 0; $i2 < $n; $i2++) {
$C[$i1][$i2] = 0.0;
for ($i3 = 0; $i3 < $n; $i3++) {
$x = ($i1 == $i3) ? $A[$i1][$i3] - $l2 : $A[$i1][$i3];
$C[$i1][$i2] += $x * $B[$i3][$i2];
}
}
}
for ($i1 = 0; $i1 < $n; $i1++) {
for ($i2 = 0; $i2 < $n; $i2++)
$B[$i1][$i2] = $C[$i1][$i2];
}
for ($i1 = 0; $i1 < $n; $i1++) {
$v1[$i1] = 0.0;
for ($i2 = 0; $i2 < $n; $i2++)
$v1[$i1] += $B[$i1][$i2] * $v0[$i2];
}
for ($i1 = 0; $i1 < $n; $i1++)
$v0[$i1] = $v1[$i1];
$ind = power($i+1, $n, $m, $ct, $eps, $A, $B, $C, $r, $v, $v0, $v1, $v2);
}
}
// 収束しない場合
else {
for ($i1 = 0; $i1 < $n; $i1++)
$v1[$i1] = $v2[$i1];
$l1 = $l2;
$k++;
}
}
return $ind;
}
/*
---------データ例（コメント部分を除いて下さい）---------
2 2 100 // 各組の変数の数（r と s， r ≦ s）とデータの数（n）
66 22 44 31 // x1, x2, x3, x4
25 74 17 81
50 23 53 71
25 57 19 81
74 47 64 47
39 33 48 46
14 22 9 69
67 60 49 26
42 40 77 65
11 80 0 86
32 0 43 74
68 69 44 68
24 49 9 71
42 74 28 46
60 58 73 28
36 37 33 68
24 44 19 83
30 40 31 50
55 40 60 49
63 47 94 41
72 30 100 45
19 22 13 75
43 39 43 34
90 83 92 31
51 77 52 82
53 70 34 31
28 51 53 44
40 62 42 79
31 48 22 68
57 29 51 30
64 89 57 42
49 82 72 29
53 31 55 43
79 52 70 10
45 19 43 57
35 34 34 89
4 69 0 100
49 49 66 66
92 82 97 6
5 89 0 100
65 26 83 28
56 36 64 38
48 50 25 22
30 30 15 55
40 65 38 42
14 67 9 67
84 96 90 8
53 64 51 54
50 89 60 52
76 41 68 9
49 40 53 49
78 66 66 17
76 58 90 29
41 15 40 49
63 60 55 33
40 36 49 67
78 54 71 18
62 72 69 12
64 47 42 53
56 64 9 15
77 35 56 25
44 12 46 87
80 9 56 19
36 21 52 78
48 63 64 48
43 61 50 47
58 23 28 50
90 12 100 0
13 33 11 77
67 44 48 28
75 45 68 17
81 22 89 9
46 45 59 55
56 49 64 55
65 62 72 27
34 49 29 77
45 33 60 63
20 45 14 99
33 38 26 87
44 51 69 52
64 57 64 48
44 64 51 28
63 48 56 11
29 39 33 84
40 48 51 54
40 38 26 62
68 46 61 26
58 45 68 48
64 44 77 63
59 62 44 66
81 53 93 19
23 34 12 68
51 35 55 46
74 70 84 17
42 33 56 44
46 31 46 53
33 57 38 63
40 24 20 42
53 36 60 31
0 34 0 100
*/
?>
############################
# 正準相関分析
# coded by Y.Suganuma
############################
############################################
# 線形連立方程式を解く（逆行列を求める）
# w : 方程式の左辺及び右辺
# n : 方程式の数
# m : 方程式の右辺の列の数
# eps : 逆行列の存在を判定する規準
# return : =0 : 正常
# =1 : 逆行列が存在しない
# coded by Y.Suganuma
############################################
def gauss(w, n, m, eps)
nm = n + m;
ind = 0
for i1 in 0 ... n
y1 = 0.0
m1 = i1 + 1
m2 = 0
# ピボット要素の選択
for i2 in i1 ... n
y2 = w[i2][i1].abs()
if y1 < y2
y1 = y2
m2 = i2
end
end
# 逆行列が存在しない
if y1 < eps
ind = 1
break
# 逆行列が存在する
else # 行の入れ替え
for i2 in i1 ... nm
y1 = w[i1][i2]
w[i1][i2] = w[m2][i2]
w[m2][i2] = y1
end
# 掃き出し操作
y1 = 1.0 / w[i1][i1]
for i2 in m1 ... nm
w[i1][i2] *= y1
end
for i2 in 0 ... n
if i2 != i1
for i3 in m1 ... nm
w[i2][i3] -= (w[i2][i1] * w[i1][i3])
end
end
end
end
end
return ind
end
#***************************************/
# 最大固有値と固有ベクトル（べき乗法） */
# i : 何番目の固有値かを示す */
# n : 次数 */
# m : 丸め誤差の修正回数 */
# ct : 最大繰り返し回数 */
# eps : 収束判定条件 */
# a : 対象とする行列 */
# b : nxnの行列，最初は，単位行列 */
# c : 作業域，nxnの行列 */
# r : 固有値 */
# v : 各行が固有ベクトル（nxn） */
# v0 : 固有ベクトルの初期値 */
# v1,v2 : 作業域（n次元ベクトル） */
# return : 求まった固有値の数 */
# coded by Y.Suganuma */
#***************************************/
def power(i, n, m, ct, eps, a, b, c, r, v, v0, v1, v2)
# 初期設定
ind = i
k = 0
l1 = 0.0
for i1 in 0 ... n
l1 += v0[i1] * v0[i1]
v1[i1] = v0[i1]
end
l1 = Math.sqrt(l1)
# 繰り返し計算
while k < ct
# 丸め誤差の修正
if k%m == 0
l2 = 0.0
for i1 in 0 ... n
v2[i1] = 0.0
for i2 in 0 ... n
v2[i1] += b[i1][i2] * v1[i2]
end
l2 += v2[i1] * v2[i1]
end
l2 = Math.sqrt(l2)
for i1 in 0 ... n
v1[i1] = v2[i1] / l2
end
end
# 次の近似
l2 = 0.0
for i1 in 0 ... n
v2[i1] = 0.0
for i2 in 0 ... n
v2[i1] += a[i1][i2] * v1[i2]
end
l2 += v2[i1] * v2[i1]
end
l2 = Math.sqrt(l2)
for i1 in 0 ... n
v2[i1] /= l2
end
# 収束判定
# 収束した場合
if ((l2-l1)/l1).abs() < eps
k1 = -1
for i1 in 0 ... n
if v2[i1].abs() > 0.001
k1 = i1
if v2[k1]*v1[k1] < 0.0
l2 = -l2
end
end
if k1 >= 0
break
end
end
k = ct
r[i] = l2
for i1 in 0 ... n
v[i][i1] = v2[i1]
end
if i == n-1
ind = i + 1
else
for i1 in 0 ... n
for i2 in 0 ... n
c[i1][i2] = 0.0
for i3 in 0 ... n
x = (i1 == i3) ? a[i1][i3] - l2 : a[i1][i3]
c[i1][i2] += x * b[i3][i2]
end
end
end
for i1 in 0 ... n
for i2 in 0 ... n
b[i1][i2] = c[i1][i2]
end
end
for i1 in 0 ... n
v1[i1] = 0.0
for i2 in 0 ... n
v1[i1] += b[i1][i2] * v0[i2]
end
end
for i1 in 0 ... n
v0[i1] = v1[i1]
end
ind = power(i+1, n, m, ct, eps, a, b, c, r, v, v0, v1, v2)
end
# 収束しない場合
else
for i1 in 0 ... n
v1[i1] = v2[i1]
end
l1 = l2
k += 1
end
end
return ind
end
###################################
# 正準相関分析
# r,s : 各組における変数の数
# n : データの数
# x : データ
# ryz : 相関係数
# ab : 各組の係数（a,bの順）
# eps : 正則性を判定する規準
# v0 : 固有ベクトルの初期値
# m : 丸め誤差の修正回数
# ct : 最大繰り返し回数
# return : >0 : 相関係数の数
# =0 : エラー
# coded by Y.Suganuma
###################################
def canonical(r, s, n, x, ryz, ab, eps, v0, m, ct)
sw = 0
# 領域の確保
q = r + s
mean = Array.new(q)
w1 = Array.new(s)
v1 = Array.new(r)
v2 = Array.new(r)
a = Array.new(s)
b = Array.new(r)
c = Array.new(q)
c11 = Array.new(r)
c11i = Array.new(r)
c12 = Array.new(r)
c21 = Array.new(s)
c22 = Array.new(s)
c22i = Array.new(s)
w = Array.new(s)
for i1 in 0 ... q
c[i1] = Array.new(q)
end
for i1 in 0 ... s
a[i1] = Array.new(s)
c21[i1] = Array.new(r)
c22[i1] = Array.new(s)
c22i[i1] = Array.new(s)
w[i1] = Array.new(2*s)
end
for i1 in 0 ... r
b[i1] = Array.new(r)
c11[i1] = Array.new(r)
c11i[i1] = Array.new(r)
c12[i1] = Array.new(s)
end
# 平均値の計算
for i1 in 0 ... q
mean[i1] = 0.0
for i2 in 0 ... n
mean[i1] += x[i1][i2]
end
mean[i1] /= n
end
# 分散強分散行列の計算
for i1 in 0 ... q
for i2 in i1 ... q
vv = 0.0
for i3 in 0 ... n
vv += (x[i1][i3] - mean[i1]) * (x[i2][i3] - mean[i2])
end
vv /= (n - 1)
c[i1][i2] = vv
if i1 != i2
c[i2][i1] = vv
end
end
end
# c11, c12, c21, c22 の設定
# c12
for i1 in 0 ... r
for i2 in 0 ... s
c12[i1][i2] = c[i1][i2+r]
end
end
# c21
for i1 in 0 ... s
for i2 in 0 ... r
c21[i1][i2] = c[i1+r][i2]
end
end
# c11とその逆行列
for i1 in 0 ... r
for i2 in 0 ... r
w[i1][i2] = c[i1][i2]
c11[i1][i2] = c[i1][i2]
end
for i2 in r ... 2*r
w[i1][i2] = 0.0
end
w[i1][i1+r] = 1.0
end
sw1 = gauss(w, r, r, 1.0e-10)
if sw1 == 0
for i1 in 0 ... r
for i2 in 0 ... r
c11i[i1][i2] = w[i1][i2+r]
end
end
else
sw = 1
end
# c22とその逆行列
for i1 in 0 ... s
for i2 in 0 ... s
w[i1][i2] = c[i1+r][i2+r]
c22[i1][i2] = c[i1+r][i2+r]
end
for i2 in s ... 2*s
w[i1][i2] = 0.0
end
w[i1][i1+s] = 1.0
end
sw1 = gauss(w, s, s, eps)
if sw1 == 0
for i1 in 0 ... s
for i2 in 0 ... s
c22i[i1][i2] = w[i1][i2+s]
end
end
else
sw = 1
end
# 固有値λ及び固有ベクトルaの計算
if sw == 0
# 行列の計算
for i1 in 0 ... s
for i2 in 0 ... r
a[i1][i2] = 0.0
for i3 in 0 ... s
a[i1][i2] += c22i[i1][i3] * c21[i3][i2]
end
end
end
for i1 in 0 ... r
for i2 in 0 ... s
w[i1][i2] = 0.0
for i3 in 0 ... s
w[i1][i2] += c12[i1][i3] * a[i3][i2]
end
end
end
for i1 in 0 ... r
for i2 in 0 ... r
a[i1][i2] = 0.0
for i3 in 0 ... r
a[i1][i2] += c11i[i1][i3] * w[i3][i2]
end
end
end
# 固有値と固有ベクトル（べき乗法）
for i1 in 0 ... r
for i2 in 0 ... r
b[i1][i2] = 0.0
end
b[i1][i1] = 1.0
end
sw = power(0, r, m, ct, eps, a, b, c, ryz, ab, v0, v1, v2)
if sw > 0
for i1 in 0 ... r
# 相関係数
ryz[i1] = Math.sqrt(ryz[i1])
# 大きさの調整(a)
for i2 in 0 ... r
w1[i2] = 0.0
for i3 in 0 ... r
w1[i2] += c11[i2][i3] * ab[i1][i3]
end
end
len = 0.0
for i2 in 0 ... r
len += ab[i1][i2] * w1[i2]
end
len = Math.sqrt(len)
for i2 in 0 ... r
ab[i1][i2] /= len
end
# bの計算
for i2 in 0 ... s
w1[i2] = 0.0
for i3 in 0 ... r
w1[i2] += c21[i2][i3] * ab[i1][i3]
end
end
for i2 in 0 ... s
ab[i1][i2+r] = 0.0
for i3 in 0 ... s
ab[i1][i2+r] += c22i[i2][i3] * w1[i3]
end
end
for i2 in 0 ... s
ab[i1][i2+r] /= ryz[i1]
end
# 大きさの調整(b)
for i2 in 0 ... s
w1[i2] = 0.0
for i3 in 0 ... s
w1[i2] += c22[i2][i3] * ab[i1][i3+r]
end
end
len = 0.0
for i2 in 0 ... s
len += ab[i1][i2+r] * w1[i2]
end
len = Math.sqrt(len)
for i2 in 0 ... s
ab[i1][i2+r] /= len
end
end
end
else
sw = 0
end
return sw
end
ss = gets().split(" ")
r = Integer(ss[0]) # 各組の変数
s = Integer(ss[1]) # 各組の変数
n = Integer(ss[2]) # データの数
q = r + s
ryz = Array.new(r)
v0 = Array.new(r)
for i1 in 0 ... r
v0[i1] = 1.0
end
x = Array.new(q)
for i1 in 0 ... q
x[i1] = Array.new(n)
end
ab = Array.new(r)
for i1 in 0 ... r
ab[i1] = Array.new(q)
end
for i1 in 0 ... n # データ
ss = gets().split()
for i2 in 0 ... q
x[i2][i1] = Float(ss[i2])
end
end
sw = canonical(r, s, n, x, ryz, ab, 1.0e-10, v0, 15, 200)
if sw > 0
for i1 in 0 ... sw
print("相関係数 " + String(ryz[i1]) + "\n")
print(" a")
for i2 in 0 ... r
print(" " + String(ab[i1][i2]))
end
print("\n")
print(" b")
for i2 in 0 ... s
print(" " + String(ab[i1][r+i2]))
end
print("\n")
end
else
print("***error 解を求めることができませんでした\n")
end
=begin
---------データ例（コメント部分を除いて下さい）---------
2 2 100 # 各組の変数の数（r と s， r ≦ s）とデータの数（n）
66 22 44 31 # x1, x2, x3, x4
25 74 17 81
50 23 53 71
25 57 19 81
74 47 64 47
39 33 48 46
14 22 9 69
67 60 49 26
42 40 77 65
11 80 0 86
32 0 43 74
68 69 44 68
24 49 9 71
42 74 28 46
60 58 73 28
36 37 33 68
24 44 19 83
30 40 31 50
55 40 60 49
63 47 94 41
72 30 100 45
19 22 13 75
43 39 43 34
90 83 92 31
51 77 52 82
53 70 34 31
28 51 53 44
40 62 42 79
31 48 22 68
57 29 51 30
64 89 57 42
49 82 72 29
53 31 55 43
79 52 70 10
45 19 43 57
35 34 34 89
4 69 0 100
49 49 66 66
92 82 97 6
5 89 0 100
65 26 83 28
56 36 64 38
48 50 25 22
30 30 15 55
40 65 38 42
14 67 9 67
84 96 90 8
53 64 51 54
50 89 60 52
76 41 68 9
49 40 53 49
78 66 66 17
76 58 90 29
41 15 40 49
63 60 55 33
40 36 49 67
78 54 71 18
62 72 69 12
64 47 42 53
56 64 9 15
77 35 56 25
44 12 46 87
80 9 56 19
36 21 52 78
48 63 64 48
43 61 50 47
58 23 28 50
90 12 100 0
13 33 11 77
67 44 48 28
75 45 68 17
81 22 89 9
46 45 59 55
56 49 64 55
65 62 72 27
34 49 29 77
45 33 60 63
20 45 14 99
33 38 26 87
44 51 69 52
64 57 64 48
44 64 51 28
63 48 56 11
29 39 33 84
40 48 51 54
40 38 26 62
68 46 61 26
58 45 68 48
64 44 77 63
59 62 44 66
81 53 93 19
23 34 12 68
51 35 55 46
74 70 84 17
42 33 56 44
46 31 46 53
33 57 38 63
40 24 20 42
53 36 60 31
0 34 0 100
=end
# -*- coding: UTF-8 -*-
import numpy as np
import sys
from math import *
############################################
# 線形連立方程式を解く（逆行列を求める）
# w : 方程式の左辺及び右辺
# n : 方程式の数
# m : 方程式の右辺の列の数
# eps : 逆行列の存在を判定する規準
# return : =0 : 正常
# =1 : 逆行列が存在しない
# coded by Y.Suganuma
############################################
def gauss(w, n, m, eps) :
nm = n + m
ind = 0
for i1 in range(0, n) :
y1 = 0.0
m1 = i1 + 1
m2 = 0
# ピボット要素の選択
for i2 in range(i1, n) :
y2 = abs(w[i2][i1])
if y1 < y2 :
y1 = y2
m2 = i2
# 逆行列が存在しない
if y1 < eps :
ind = 1
break
# 逆行列が存在する
else : # 行の入れ替え
for i2 in range(i1, nm) :
y1 = w[i1][i2]
w[i1][i2] = w[m2][i2]
w[m2][i2] = y1
# 掃き出し操作
y1 = 1.0 / w[i1][i1]
for i2 in range(m1, nm) :
w[i1][i2] *= y1
for i2 in range(0, n) :
if i2 != i1 :
for i3 in range(m1, nm) :
w[i2][i3] -= (w[i2][i1] * w[i1][i3])
return ind
############################################
# 最大固有値と固有ベクトル（べき乗法）
# i : 何番目の固有値かを示す
# n : 次数
# m : 丸め誤差の修正回数
# ct : 最大繰り返し回数
# eps : 収束判定条件
# A : 対象とする行列
# B : nxnの行列，最初は，単位行列
# C : 作業域，nxnの行列
# r : 固有値
# v : 各行が固有ベクトル（nxn）
# v0 : 固有ベクトルの初期値
# v1,v2 : 作業域（n次元ベクトル）
# return : 求まった固有値の数
# coded by Y.Suganuma
############################################
def power(i, n, m, ct, eps, A, B, C, r, v, v0, v1, v2) :
# 初期設定
ind = i
k = 0
l1 = 0.0
for i1 in range(0, n) :
l1 += v0[i1] * v0[i1]
v1[i1] = v0[i1]
l1 = sqrt(l1)
# 繰り返し計算
while k < ct :
# 丸め誤差の修正
if k%m == 0 :
l2 = 0.0
for i1 in range(0, n) :
v2[i1] = 0.0
for i2 in range(0, n) :
v2[i1] += B[i1][i2] * v1[i2]
l2 += v2[i1] * v2[i1]
l2 = sqrt(l2)
for i1 in range(0, n) :
v1[i1] = v2[i1] / l2
# 次の近似
l2 = 0.0
for i1 in range(0, n) :
v2[i1] = 0.0
for i2 in range(0, n) :
v2[i1] += A[i1][i2] * v1[i2]
l2 += v2[i1] * v2[i1]
l2 = sqrt(l2)
for i1 in range(0, n) :
v2[i1] /= l2
# 収束判定
# 収束した場合
if abs((l2-l1)/l1) < eps :
k1 = -1
for i1 in range(0, n) :
if abs(v2[i1]) > 0.001 :
k1 = i1
if v2[k1]*v1[k1] < 0.0 :
l2 = -l2
break
k = ct
r[i] = l2
for i1 in range(0, n) :
v[i][i1] = v2[i1]
if i == n-1 :
ind = i + 1
else :
for i1 in range(0, n) :
for i2 in range(0, n) :
C[i1][i2] = 0.0
for i3 in range(0, n) :
if i1 == i3 :
x = A[i1][i3] - l2
else :
x = A[i1][i3]
C[i1][i2] += x * B[i3][i2]
for i1 in range(0, n) :
for i2 in range(0, n) :
B[i1][i2] = C[i1][i2]
for i1 in range(0, n) :
v1[i1] = 0.0
for i2 in range(0, n) :
v1[i1] += B[i1][i2] * v0[i2]
for i1 in range(0, n) :
v0[i1] = v1[i1]
ind = power(i+1, n, m, ct, eps, A, B, C, r, v, v0, v1, v2)
# 収束しない場合
else :
for i1 in range(0, n) :
v1[i1] = v2[i1]
l1 = l2
k += 1
return ind
###################################
# 正準相関分析
# r,s : 各組における変数の数
# n : データの数
# x : データ
# ryz : 相関係数
# ab : 各組の係数（a,bの順）
# eps : 正則性を判定する規準
# v0 : 固有ベクトルの初期値
# m : 丸め誤差の修正回数
# ct : 最大繰り返し回数
# return : >0 : 相関係数の数
# =0 : エラー
# coded by Y.Suganuma
###################################
def canonical(r, s, n, x, ryz, ab, eps, v0, m, ct) :
sw = 0
# 領域の確保
q = r + s
mean = np.empty(q, np.float)
w1 = np.empty(s, np.float)
v1 = np.empty(r, np.float)
v2 = np.empty(r, np.float)
A = np.empty((s, s), np.float)
B = np.empty((r, r), np.float)
C = np.empty((q, q), np.float)
C11 = np.empty((r, r), np.float)
C11i = np.empty((r, r), np.float)
C12 = np.empty((r, s), np.float)
C21 = np.empty((s, r), np.float)
C22 = np.empty((s, s), np.float)
C22i = np.empty((s, s), np.float)
w = np.empty((s, 2*s), np.float)
# 平均値の計算
for i1 in range(0, q) :
mean[i1] = 0.0
for i2 in range(0, n) :
mean[i1] += x[i1][i2]
mean[i1] /= n
# 分散強分散行列の計算
for i1 in range(0, q) :
for i2 in range(i1, q) :
vv = 0.0
for i3 in range(0, n) :
vv += (x[i1][i3] - mean[i1]) * (x[i2][i3] - mean[i2])
vv /= (n - 1)
C[i1][i2] = vv
if i1 != i2 :
C[i2][i1] = vv
# C11, C12, C21, C22 の設定
# C12
for i1 in range(0, r) :
for i2 in range(0, s) :
C12[i1][i2] = C[i1][i2+r]
# C21
for i1 in range(0, s) :
for i2 in range(0, r) :
C21[i1][i2] = C[i1+r][i2]
# C11とその逆行列
for i1 in range(0, r) :
for i2 in range(0, r) :
w[i1][i2] = C[i1][i2]
C11[i1][i2] = C[i1][i2]
for i2 in range(r, 2*r) :
w[i1][i2] = 0.0
w[i1][i1+r] = 1.0
sw1 = gauss(w, r, r, 1.0e-10)
if sw1 == 0 :
for i1 in range(0, r) :
for i2 in range(0, r) :
C11i[i1][i2] = w[i1][i2+r]
else :
sw = 1
# C22とその逆行列
for i1 in range(0, s) :
for i2 in range(0, s) :
w[i1][i2] = C[i1+r][i2+r]
C22[i1][i2] = C[i1+r][i2+r]
for i2 in range(s, 2*s) :
w[i1][i2] = 0.0
w[i1][i1+s] = 1.0
sw1 = gauss(w, s, s, eps)
if sw1 == 0 :
for i1 in range(0, s) :
for i2 in range(0, s) :
C22i[i1][i2] = w[i1][i2+s]
else :
sw = 1
# 固有値λ及び固有ベクトルaの計算
if sw == 0 :
# 行列の計算
for i1 in range(0, s) :
for i2 in range(0, r) :
A[i1][i2] = 0.0
for i3 in range(0, s) :
A[i1][i2] += C22i[i1][i3] * C21[i3][i2]
for i1 in range(0, r) :
for i2 in range(0, s) :
w[i1][i2] = 0.0
for i3 in range(0, s) :
w[i1][i2] += C12[i1][i3] * A[i3][i2]
for i1 in range(0, r) :
for i2 in range(0, r) :
A[i1][i2] = 0.0
for i3 in range(0, r) :
A[i1][i2] += C11i[i1][i3] * w[i3][i2]
# 固有値と固有ベクトル（べき乗法）
for i1 in range(0, r) :
for i2 in range(0, r) :
B[i1][i2] = 0.0
B[i1][i1] = 1.0
sw = power(0, r, m, ct, eps, A, B, C, ryz, ab, v0, v1, v2)
if sw > 0 :
for i1 in range(0, r) :
# 相関係数
ryz[i1] = sqrt(ryz[i1])
# 大きさの調整(a)
for i2 in range(0, r) :
w1[i2] = 0.0
for i3 in range(0, r) :
w1[i2] += C11[i2][i3] * ab[i1][i3]
len = 0.0
for i2 in range(0, r) :
len += ab[i1][i2] * w1[i2]
len = sqrt(len)
for i2 in range(0, r) :
ab[i1][i2] /= len
# bの計算
for i2 in range(0, s) :
w1[i2] = 0.0
for i3 in range(0, r) :
w1[i2] += C21[i2][i3] * ab[i1][i3]
for i2 in range(0, s) :
ab[i1][i2+r] = 0.0
for i3 in range(0, s) :
ab[i1][i2+r] += C22i[i2][i3] * w1[i3]
for i2 in range(0, s) :
ab[i1][i2+r] /= ryz[i1]
# 大きさの調整(b)
for i2 in range(0, s) :
w1[i2] = 0.0
for i3 in range(0, s) :
w1[i2] += C22[i2][i3] * ab[i1][i3+r]
len = 0.0
for i2 in range(0, s) :
len += ab[i1][i2+r] * w1[i2]
len = sqrt(len)
for i2 in range(0, s) :
ab[i1][i2+r] /= len
else :
sw = 0
return sw
############################
# 正準相関分析
# coded by Y.Suganuma
############################
line = sys.stdin.readline()
ss = line.split()
r = int(ss[0]) # 各組の変数
s = int(ss[1]) # 各組の変数
n = int(ss[2]) # データの数
q = r + s
ryz = np.empty(r, np.float)
v0 = np.ones(r, np.float)
x = np.empty((q, n), np.float)
ab = np.empty((r, q), np.float)
for i1 in range(0, n) : # データ
line = sys.stdin.readline()
ss = line.split()
for i2 in range(0, q) :
x[i2][i1] = float(ss[i2])
sw = canonical(r, s, n, x, ryz, ab, 1.0e-10, v0, 15, 200)
if sw > 0 :
for i1 in range(0, sw) :
print("相関係数 " + str(ryz[i1]))
print(" a", end="")
for i2 in range(0, r) :
print(" " + str(ab[i1][i2]), end="")
print()
print(" b", end="")
for i2 in range(0, s) :
print(" " + str(ab[i1][r+i2]), end="")
print()
else :
print("***error 解を求めることができませんでした")
"""
---------データ例（コメント部分を除いて下さい）---------
2 2 100 # 各組の変数の数（r と s， r ≦ s）とデータの数（n）
66 22 44 31 # x1, x2, x3, x4
25 74 17 81
50 23 53 71
25 57 19 81
74 47 64 47
39 33 48 46
14 22 9 69
67 60 49 26
42 40 77 65
11 80 0 86
32 0 43 74
68 69 44 68
24 49 9 71
42 74 28 46
60 58 73 28
36 37 33 68
24 44 19 83
30 40 31 50
55 40 60 49
63 47 94 41
72 30 100 45
19 22 13 75
43 39 43 34
90 83 92 31
51 77 52 82
53 70 34 31
28 51 53 44
40 62 42 79
31 48 22 68
57 29 51 30
64 89 57 42
49 82 72 29
53 31 55 43
79 52 70 10
45 19 43 57
35 34 34 89
4 69 0 100
49 49 66 66
92 82 97 6
5 89 0 100
65 26 83 28
56 36 64 38
48 50 25 22
30 30 15 55
40 65 38 42
14 67 9 67
84 96 90 8
53 64 51 54
50 89 60 52
76 41 68 9
49 40 53 49
78 66 66 17
76 58 90 29
41 15 40 49
63 60 55 33
40 36 49 67
78 54 71 18
62 72 69 12
64 47 42 53
56 64 9 15
77 35 56 25
44 12 46 87
80 9 56 19
36 21 52 78
48 63 64 48
43 61 50 47
58 23 28 50
90 12 100 0
13 33 11 77
67 44 48 28
75 45 68 17
81 22 89 9
46 45 59 55
56 49 64 55
65 62 72 27
34 49 29 77
45 33 60 63
20 45 14 99
33 38 26 87
44 51 69 52
64 57 64 48
44 64 51 28
63 48 56 11
29 39 33 84
40 48 51 54
40 38 26 62
68 46 61 26
58 45 68 48
64 44 77 63
59 62 44 66
81 53 93 19
23 34 12 68
51 35 55 46
74 70 84 17
42 33 56 44
46 31 46 53
33 57 38 63
40 24 20 42
53 36 60 31
0 34 0 100
"""
/****************************/
/* 正準相関分析 */
/* coded by Y.Suganuma */
/****************************/
using System;
class Program
{
static void Main()
{
Test1 ts = new Test1();
}
}
class Test1
{
public Test1()
{
char[] charSep = new char[] {' '};
// 各組の変数の数とデータの数
String[] str = Console.ReadLine().Split(charSep, StringSplitOptions.RemoveEmptyEntries);
int r = int.Parse(str[0]);
int s = int.Parse(str[1]);
int n = int.Parse(str[2]);
int q = r + s;
double[] ryz = new double [r];
double[] v0 = new double [r];
for (int i1 = 0; i1 < r; i1++)
v0[i1] = 1.0;
double[][] x = new double [q][];
for (int i1 = 0; i1 < q; i1++)
x[i1] = new double [n];
double[][] ab = new double [r][];
for (int i1 = 0; i1 < r; i1++)
ab[i1] = new double [q];
// データ
for (int i1 = 0; i1 < n; i1++) {
str = Console.ReadLine().Split(charSep, StringSplitOptions.RemoveEmptyEntries);
for (int i2 = 0; i2 < q; i2++)
x[i2][i1] = int.Parse(str[i2]);
}
int sw = canonical(r, s, n, x, ryz, ab, 1.0e-10, v0, 15, 200);
if (sw > 0) {
for (int i1 = 0; i1 < sw; i1++) {
Console.WriteLine("相関係数 " + ryz[i1]);
Console.Write(" a");
for (int i2 = 0; i2 < r; i2++)
Console.Write(" " + ab[i1][i2]);
Console.WriteLine();
Console.Write(" b");
for (int i2 = 0; i2 < s; i2++)
Console.Write(" " + ab[i1][r+i2]);
Console.WriteLine();
}
}
else
Console.WriteLine("***error 解を求めることができませんでした");
}
/***********************************/
/* 正準相関分析 */
/* r,s : 各組における変数の数 */
/* n : データの数 */
/* x : データ */
/* ryz : 相関係数 */
/* ab : 各組の係数（a,bの順） */
/* eps : 正則性を判定する規準 */
/* v0 : 固有ベクトルの初期値 */
/* m : 丸め誤差の修正回数 */
/* ct : 最大繰り返し回数 */
/* return : >0 : 相関係数の数 */
/* =0 : エラー */
/***********************************/
int canonical(int r, int s, int n, double[][] x, double[] ryz,
double[][] ab, double eps, double[] v0, int m, int ct)
{
int sw = 0;
// 領域の確保
int q = r + s;
double[] mean = new double [q];
double[] w1 = new double [s];
double[] v1 = new double [r];
double[] v2 = new double [r];
double[][] A = new double [s][];
double[][] B = new double [r][];
double[][] C = new double [q][];
double[][] C11 = new double [r][];
double[][] C11i = new double [r][];
double[][] C12 = new double [r][];
double[][] C21 = new double [s][];
double[][] C22 = new double [s][];
double[][] C22i = new double [s][];
double[][] w = new double [s][];
for (int i1 = 0; i1 < q; i1++)
C[i1] = new double [q];
for (int i1 = 0; i1 < s; i1++) {
A[i1] = new double [s];
C21[i1] = new double [r];
C22[i1] = new double [s];
C22i[i1] = new double [s];
w[i1] = new double [2*s];
}
for (int i1 = 0; i1 < r; i1++) {
B[i1] = new double [r];
C11[i1] = new double [r];
C11i[i1] = new double [r];
C12[i1] = new double [s];
}
// 平均値の計算
for (int i1 = 0; i1 < q; i1++) {
mean[i1] = 0.0;
for (int i2 = 0; i2 < n; i2++)
mean[i1] += x[i1][i2];
mean[i1] /= n;
}
// 分散強分散行列の計算
for (int i1 = 0; i1 < q; i1++) {
for (int i2 = i1; i2 < q; i2++) {
double vv = 0.0;
for (int i3 = 0; i3 < n; i3++)
vv += (x[i1][i3] - mean[i1]) * (x[i2][i3] - mean[i2]);
vv /= (n - 1);
C[i1][i2] = vv;
if (i1 != i2)
C[i2][i1] = vv;
}
}
// C11, C12, C21, C22 の設定
// C12
for (int i1 = 0; i1 < r; i1++) {
for (int i2 = 0; i2 < s; i2++)
C12[i1][i2] = C[i1][i2+r];
}
// C21
for (int i1 = 0; i1 < s; i1++) {
for (int i2 = 0; i2 < r; i2++)
C21[i1][i2] = C[i1+r][i2];
}
// C11とその逆行列
for (int i1 = 0; i1 < r; i1++) {
for (int i2 = 0; i2 < r; i2++) {
w[i1][i2] = C[i1][i2];
C11[i1][i2] = C[i1][i2];
}
for (int i2 = r; i2 < 2*r; i2++)
w[i1][i2] = 0.0;
w[i1][i1+r] = 1.0;
}
int sw1 = gauss(w, r, r, 1.0e-10);
if (sw1 == 0) {
for (int i1 = 0; i1 < r; i1++) {
for (int i2 = 0; i2 < r; i2++)
C11i[i1][i2] = w[i1][i2+r];
}
}
else
sw = 1;
// C22とその逆行列
for (int i1 = 0; i1 < s; i1++) {
for (int i2 = 0; i2 < s; i2++) {
w[i1][i2] = C[i1+r][i2+r];
C22[i1][i2] = C[i1+r][i2+r];
}
for (int i2 = s; i2 < 2*s; i2++)
w[i1][i2] = 0.0;
w[i1][i1+s] = 1.0;
}
sw1 = gauss(w, s, s, eps);
if (sw1 == 0) {
for (int i1 = 0; i1 < s; i1++) {
for (int i2 = 0; i2 < s; i2++)
C22i[i1][i2] = w[i1][i2+s];
}
}
else
sw = 1;
// 固有値λ及び固有ベクトルaの計算
if (sw == 0) {
// 行列の計算
for (int i1 = 0; i1 < s; i1++) {
for (int i2 = 0; i2 < r; i2++) {
A[i1][i2] = 0.0;
for (int i3 = 0; i3 < s; i3++)
A[i1][i2] += C22i[i1][i3] * C21[i3][i2];
}
}
for (int i1 = 0; i1 < r; i1++) {
for (int i2 = 0; i2 < s; i2++) {
w[i1][i2] = 0.0;
for (int i3 = 0; i3 < s; i3++)
w[i1][i2] += C12[i1][i3] * A[i3][i2];
}
}
for (int i1 = 0; i1 < r; i1++) {
for (int i2 = 0; i2 < r; i2++) {
A[i1][i2] = 0.0;
for (int i3 = 0; i3 < r; i3++)
A[i1][i2] += C11i[i1][i3] * w[i3][i2];
}
}
// 固有値と固有ベクトル（べき乗法）
for (int i1 = 0; i1 < r; i1++) {
for (int i2 = 0; i2 < r; i2++)
B[i1][i2] = 0.0;
B[i1][i1] = 1.0;
}
sw = power(0, r, m, ct, eps, A, B, C, ryz, ab, v0, v1, v2);
if (sw > 0) {
for (int i1 = 0; i1 < r; i1++) {
// 相関係数
ryz[i1] = Math.Sqrt(ryz[i1]);
// 大きさの調整(a)
for (int i2 = 0; i2 < r; i2++) {
w1[i2] = 0.0;
for (int i3 = 0; i3 < r; i3++)
w1[i2] += C11[i2][i3] * ab[i1][i3];
}
double len = 0.0;
for (int i2 = 0; i2 < r; i2++)
len += ab[i1][i2] * w1[i2];
len = Math.Sqrt(len);
for (int i2 = 0; i2 < r; i2++)
ab[i1][i2] /= len;
// bの計算
for (int i2 = 0; i2 < s; i2++) {
w1[i2] = 0.0;
for (int i3 = 0; i3 < r; i3++)
w1[i2] += C21[i2][i3] * ab[i1][i3];
}
for (int i2 = 0; i2 < s; i2++) {
ab[i1][i2+r] = 0.0;
for (int i3 = 0; i3 < s; i3++)
ab[i1][i2+r] += C22i[i2][i3] * w1[i3];
}
for (int i2 = 0; i2 < s; i2++)
ab[i1][i2+r] /= ryz[i1];
// 大きさの調整(b)
for (int i2 = 0; i2 < s; i2++) {
w1[i2] = 0.0;
for (int i3 = 0; i3 < s; i3++)
w1[i2] += C22[i2][i3] * ab[i1][i3+r];
}
len = 0.0;
for (int i2 = 0; i2 < s; i2++)
len += ab[i1][i2+r] * w1[i2];
len = Math.Sqrt(len);
for (int i2 = 0; i2 < s; i2++)
ab[i1][i2+r] /= len;
}
}
}
else
sw = 0;
return sw;
}
/****************************************/
/* 最大固有値と固有ベクトル（べき乗法） */
/* i : 何番目の固有値かを示す */
/* n : 次数 */
/* m : 丸め誤差の修正回数 */
/* ct : 最大繰り返し回数 */
/* eps : 収束判定条件 */
/* A : 対象とする行列 */
/* B : nxnの行列，最初は，単位行列 */
/* C : 作業域，nxnの行列 */
/* r : 固有値 */
/* v : 各行が固有ベクトル（nxn） */
/* v0 : 固有ベクトルの初期値 */
/* v1,v2 : 作業域（n次元ベクトル） */
/* return : 求まった固有値の数 */
/* coded by Y.Suganuma */
/****************************************/
int power(int i, int n, int m, int ct, double eps, double[][] A,
double[][] B, double[][] C, double[] r, double[][] v,
double[] v0, double[] v1, double[] v2)
{
// 初期設定
int ind = i;
int k = 0;
double l1 = 0.0;
for (int i1 = 0; i1 < n; i1++) {
l1 += v0[i1] * v0[i1];
v1[i1] = v0[i1];
}
l1 = Math.Sqrt(l1);
// 繰り返し計算
while (k < ct) {
// 丸め誤差の修正
double l2;
if (k%m == 0) {
l2 = 0.0;
for (int i1 = 0; i1 < n; i1++) {
v2[i1] = 0.0;
for (int i2 = 0; i2 < n; i2++)
v2[i1] += B[i1][i2] * v1[i2];
l2 += v2[i1] * v2[i1];
}
l2 = Math.Sqrt(l2);
for (int i1 = 0; i1 < n; i1++)
v1[i1] = v2[i1] / l2;
}
// 次の近似
l2 = 0.0;
for (int i1 = 0; i1 < n; i1++) {
v2[i1] = 0.0;
for (int i2 = 0; i2 < n; i2++)
v2[i1] += A[i1][i2] * v1[i2];
l2 += v2[i1] * v2[i1];
}
l2 = Math.Sqrt(l2);
for (int i1 = 0; i1 < n; i1++)
v2[i1] /= l2;
// 収束判定
// 収束した場合
if (Math.Abs((l2-l1)/l1) < eps) {
int k1 = -1;
for (int i1 = 0; i1 < n && k1 < 0; i1++) {
if (Math.Abs(v2[i1]) > 0.001) {
k1 = i1;
if (v2[k1]*v1[k1] < 0.0)
l2 = -l2;
}
}
k = ct;
r[i] = l2;
for (int i1 = 0; i1 < n; i1++)
v[i][i1] = v2[i1];
if (i == n-1)
ind = i + 1;
else {
for (int i1 = 0; i1 < n; i1++) {
for (int i2 = 0; i2 < n; i2++) {
C[i1][i2] = 0.0;
for (int i3 = 0; i3 < n; i3++) {
double x = (i1 == i3) ? A[i1][i3] - l2 : A[i1][i3];
C[i1][i2] += x * B[i3][i2];
}
}
}
for (int i1 = 0; i1 < n; i1++) {
for (int i2 = 0; i2 < n; i2++)
B[i1][i2] = C[i1][i2];
}
for (int i1 = 0; i1 < n; i1++) {
v1[i1] = 0.0;
for (int i2 = 0; i2 < n; i2++)
v1[i1] += B[i1][i2] * v0[i2];
}
for (int i1 = 0; i1 < n; i1++)
v0[i1] = v1[i1];
ind = power(i+1, n, m, ct, eps, A, B, C, r, v, v0, v1, v2);
}
}
// 収束しない場合
else {
for (int i1 = 0; i1 < n; i1++)
v1[i1] = v2[i1];
l1 = l2;
k++;
}
}
return ind;
}
/*******************************************************/
/* 線形連立方程式を解く（逆行列を求める） */
/* w : 方程式の左辺及び右辺 */
/* n : 方程式の数 */
/* m : 方程式の右辺の列の数 */
/* eps : 逆行列の存在を判定する規準 */
/* return : =0 : 正常 */
/* =1 : 逆行列が存在しない */
/*******************************************************/
int gauss(double[][] w, int n, int m, double eps)
{
int ind = 0;
int nm = n + m;
for (int i1 = 0; i1 < n && ind == 0; i1++) {
double y1 = .0;
int m1 = i1 + 1;
int m2 = 0;
// ピボット要素の選択
for (int i2 = i1; i2 < n; i2++) {
double y2 = Math.Abs(w[i2][i1]);
if (y1 < y2) {
y1 = y2;
m2 = i2;
}
}
// 逆行列が存在しない
if (y1 < eps)
ind = 1;
// 逆行列が存在する
else {
// 行の入れ替え
for (int i2 = i1; i2 < nm; i2++) {
y1 = w[i1][i2];
w[i1][i2] = w[m2][i2];
w[m2][i2] = y1;
}
// 掃き出し操作
y1 = 1.0 / w[i1][i1];
for (int i2 = m1; i2 < nm; i2++)
w[i1][i2] *= y1;
for (int i2 = 0; i2 < n; i2++) {
if (i2 != i1) {
for (int i3 = m1; i3 < nm; i3++)
w[i2][i3] -= w[i2][i1] * w[i1][i3];
}
}
}
}
return ind;
}
}
/*
---------データ例（コメント部分を除いて下さい）---------
2 2 100 // 各組の変数の数（r と s， r ≦ s）とデータの数（n）
66 22 44 31 // x1, x2, x3, x4
25 74 17 81
50 23 53 71
25 57 19 81
74 47 64 47
39 33 48 46
14 22 9 69
67 60 49 26
42 40 77 65
11 80 0 86
32 0 43 74
68 69 44 68
24 49 9 71
42 74 28 46
60 58 73 28
36 37 33 68
24 44 19 83
30 40 31 50
55 40 60 49
63 47 94 41
72 30 100 45
19 22 13 75
43 39 43 34
90 83 92 31
51 77 52 82
53 70 34 31
28 51 53 44
40 62 42 79
31 48 22 68
57 29 51 30
64 89 57 42
49 82 72 29
53 31 55 43
79 52 70 10
45 19 43 57
35 34 34 89
4 69 0 100
49 49 66 66
92 82 97 6
5 89 0 100
65 26 83 28
56 36 64 38
48 50 25 22
30 30 15 55
40 65 38 42
14 67 9 67
84 96 90 8
53 64 51 54
50 89 60 52
76 41 68 9
49 40 53 49
78 66 66 17
76 58 90 29
41 15 40 49
63 60 55 33
40 36 49 67
78 54 71 18
62 72 69 12
64 47 42 53
56 64 9 15
77 35 56 25
44 12 46 87
80 9 56 19
36 21 52 78
48 63 64 48
43 61 50 47
58 23 28 50
90 12 100 0
13 33 11 77
67 44 48 28
75 45 68 17
81 22 89 9
46 45 59 55
56 49 64 55
65 62 72 27
34 49 29 77
45 33 60 63
20 45 14 99
33 38 26 87
44 51 69 52
64 57 64 48
44 64 51 28
63 48 56 11
29 39 33 84
40 48 51 54
40 38 26 62
68 46 61 26
58 45 68 48
64 44 77 63
59 62 44 66
81 53 93 19
23 34 12 68
51 35 55 46
74 70 84 17
42 33 56 44
46 31 46 53
33 57 38 63
40 24 20 42
53 36 60 31
0 34 0 100
*/
'**************************'
' 正準相関分析 '
' coded by Y.Suganuma '
'**************************'
Imports System.Text.RegularExpressions
Module Test
Sub Main()
Dim MS As Regex = New Regex("\s+")
' 各組の変数の数とデータの数
Dim str() As String = MS.Split(Console.ReadLine().Trim())
Dim r As Integer = Integer.Parse(str(0))
Dim s As Integer = Integer.Parse(str(1))
Dim n As Integer = Integer.Parse(str(2))
Dim q As Integer = r + s
Dim ryz(r) As Double
Dim v0(r) As Double
For i1 As Integer = 0 To r-1
v0(i1) = 1.0
Next
Dim x(q,n) As Double
Dim ab(r,q) As Double
' データ
For i1 As Integer = 0 To n-1
str = MS.Split(Console.ReadLine().Trim())
For i2 As Integer = 0 To q-1
x(i2,i1) = Integer.Parse(str(i2))
Next
Next
Dim sw As Integer = canonical(r, s, n, x, ryz, ab, 1.0e-10, v0, 15, 200)
If sw > 0
For i1 As Integer = 0 To sw-1
Console.WriteLine("相関係数 " & ryz(i1))
Console.Write(" a")
For i2 As Integer = 0 To r-1
Console.Write(" " & ab(i1,i2))
Next
Console.WriteLine()
Console.Write(" b")
For i2 As Integer = 0 To s-1
Console.Write(" " & ab(i1,r+i2))
Next
Console.WriteLine()
Next
Else
Console.WriteLine("***error 解を求めることができませんでした")
End If
End Sub
'*********************************'
' 正準相関分析 '
' r,s : 各組における変数の数 '
' n : データの数 '
' x : データ '
' ryz : 相関係数 '
' ab : 各組の係数（a,bの順） '
' eps : 正則性を判定する規準 '
' v0 : 固有ベクトルの初期値 '
' m : 丸め誤差の修正回数 '
' ct : 最大繰り返し回数 '
' return : >0 : 相関係数の数 '
' =0 : エラー '
'*********************************'
Function canonical(r As Integer, s As Integer, n As Integer, x(,) As Double,
ryz() As Double, ab(,) As Double, eps As Double, v0() As Double,
m As Integer, ct As Integer)
Dim sw As Integer = 0
Dim q As Integer = r + s
' 領域の確保
Dim mean(q) As Double
Dim w1(s) As Double
Dim v1(r) As Double
Dim v2(r) As Double
Dim A(s,s) As Double
Dim B(r,r) As Double
Dim C(q,q) As Double
Dim C11(r,r) As Double
Dim C11i(r,r) As Double
Dim C12(r,s) As Double
Dim C21(s,r) As Double
Dim C22(s,s) As Double
Dim C22i(s,s) As Double
Dim w(s,2*s) As Double
' 平均値の計算
For i1 As Integer = 0 To q-1
mean(i1) = 0.0
For i2 As Integer = 0 To n-1
mean(i1) += x(i1,i2)
Next
mean(i1) /= n
Next
' 分散強分散行列の計算
For i1 As Integer = 0 To q-1
For i2 As Integer = i1 To q-1
Dim vv As Double = 0.0
For i3 As Integer = 0 To n-1
vv += (x(i1,i3) - mean(i1)) * (x(i2,i3) - mean(i2))
Next
vv /= (n - 1)
C(i1,i2) = vv
If i1 <> i2
C(i2,i1) = vv
End If
Next
Next
' C11, C12, C21, C22 の設定
' C12
For i1 As Integer = 0 To r-1
For i2 As Integer = 0 To s-1
C12(i1,i2) = C(i1,i2+r)
Next
Next
' C21
For i1 As Integer = 0 To s-1
For i2 As Integer = 0 To r-1
C21(i1,i2) = C(i1+r,i2)
Next
Next
' C11とその逆行列
For i1 As Integer = 0 To r-1
For i2 As Integer = 0 To r-1
w(i1,i2) = C(i1,i2)
C11(i1,i2) = C(i1,i2)
Next
For i2 As Integer = r To 2*r-1
w(i1,i2) = 0.0
Next
w(i1,i1+r) = 1.0
Next
Dim sw1 As Integer = gauss(w, r, r, 1.0e-10)
If sw1 = 0
For i1 As Integer = 0 To r-1
For i2 As Integer = 0 To r-1
C11i(i1,i2) = w(i1,i2+r)
Next
Next
Else
sw = 1
End If
' C22とその逆行列
For i1 As Integer = 0 To s-1
For i2 As Integer = 0 To s-1
w(i1,i2) = C(i1+r,i2+r)
C22(i1,i2) = C(i1+r,i2+r)
Next
For i2 As Integer = s To 2*s-1
w(i1,i2) = 0.0
Next
w(i1,i1+s) = 1.0
Next
sw1 = gauss(w, s, s, eps)
If sw1 = 0
For i1 As Integer = 0 To s-1
For i2 As Integer = 0 To s-1
C22i(i1,i2) = w(i1,i2+s)
Next
Next
Else
sw = 1
End If
' 固有値λ及び固有ベクトルaの計算
If sw = 0
' 行列の計算
For i1 As Integer = 0 To s-1
For i2 As Integer = 0 To r-1
A(i1,i2) = 0.0
For i3 As Integer = 0 To s-1
A(i1,i2) += C22i(i1,i3) * C21(i3,i2)
Next
Next
Next
For i1 As Integer = 0 To r-1
For i2 As Integer = 0 To s-1
w(i1,i2) = 0.0
For i3 As Integer = 0 To s-1
w(i1,i2) += C12(i1,i3) * A(i3,i2)
Next
Next
Next
For i1 As Integer = 0 To r-1
For i2 As Integer = 0 To r-1
A(i1,i2) = 0.0
For i3 As Integer = 0 To r-1
A(i1,i2) += C11i(i1,i3) * w(i3,i2)
Next
Next
Next
' 固有値と固有ベクトル（べき乗法）
For i1 As Integer = 0 To r-1
For i2 As Integer = 0 To r-1
B(i1,i2) = 0.0
Next
B(i1,i1) = 1.0
Next
sw = power(0, r, m, ct, eps, A, B, C, ryz, ab, v0, v1, v2)
If sw > 0
For i1 As Integer = 0 To r-1
' 相関係数
ryz(i1) = Math.Sqrt(ryz(i1))
' 大きさの調整(a)
For i2 As Integer = 0 To r-1
w1(i2) = 0.0
For i3 As Integer = 0 To r-1
w1(i2) += C11(i2,i3) * ab(i1,i3)
Next
Next
Dim len As Double = 0.0
For i2 As Integer = 0 To r-1
len += ab(i1,i2) * w1(i2)
Next
len = Math.Sqrt(len)
For i2 As Integer = 0 To r-1
ab(i1,i2) /= len
Next
' bの計算
For i2 As Integer = 0 To s-1
w1(i2) = 0.0
For i3 As Integer = 0 To r-1
w1(i2) += C21(i2,i3) * ab(i1,i3)
Next
Next
For i2 As Integer = 0 To s-1
ab(i1,i2+r) = 0.0
For i3 As Integer = 0To s-1
ab(i1,i2+r) += C22i(i2,i3) * w1(i3)
Next
Next
For i2 As Integer = 0 To s-1
ab(i1,i2+r) /= ryz(i1)
Next
' 大きさの調整(b)
For i2 As Integer = 0 To s-1
w1(i2) = 0.0
For i3 As Integer = 0 To s-1
w1(i2) += C22(i2,i3) * ab(i1,i3+r)
Next
Next
len = 0.0
For i2 As Integer = 0 To s-1
len += ab(i1,i2+r) * w1(i2)
Next
len = Math.Sqrt(len)
For i2 As Integer = 0 To s-1
ab(i1,i2+r) /= len
Next
Next
End If
Else
sw = 0
End If
Return sw
End Function
''''''''''''''''''''''''''''''''''''''''
' 最大固有値と固有ベクトル（べき乗法） '
' i : 何番目の固有値かを示す '
' n : 次数 '
' m : 丸め誤差の修正回数 '
' ct : 最大繰り返し回数 '
' eps : 収束判定条件 '
' A : 対象とする行列 '
' B : nxnの行列，最初は，単位行列 '
' C : 作業域，nxnの行列 '
' r : 固有値 '
' v : 各行が固有ベクトル（nxn） '
' v0 : 固有ベクトルの初期値 '
' v1,v2 : 作業域（n次元ベクトル） '
' return : 求まった固有値の数 '
' coded by Y.Suganuma '
''''''''''''''''''''''''''''''''''''''''
Function power(i As Integer, n As Integer, m As Integer, ct As Integer,
eps As Double, A(,) As Double, B(,) As Double, C(,) As Double,
r() As Double, v(,) As Double, v0() As Double, v1() As Double,
v2() As Double)
' 初期設定
Dim ind As Integer = i
Dim k As Integer = 0
Dim l1 As Double = 0.0
For i1 As Integer = 0 To n-1
l1 += v0(i1) * v0(i1)
v1(i1) = v0(i1)
Next
l1 = Math.Sqrt(l1)
' 繰り返し計算
Do While k < ct
' 丸め誤差の修正
Dim l2 As Double
If (k Mod m) = 0
l2 = 0.0
For i1 As Integer = 0 To n-1
v2(i1) = 0.0
For i2 As Integer = 0 To n-1
v2(i1) += B(i1,i2) * v1(i2)
Next
l2 += v2(i1) * v2(i1)
Next
l2 = Math.Sqrt(l2)
For i1 As Integer = 0 To n-1
v1(i1) = v2(i1) / l2
Next
End If
' 次の近似
l2 = 0.0
For i1 As Integer = 0 To n-1
v2(i1) = 0.0
For i2 As Integer = 0 To n-1
v2(i1) += A(i1,i2) * v1(i2)
Next
l2 += v2(i1) * v2(i1)
Next
l2 = Math.Sqrt(l2)
For i1 As Integer = 0 To n-1
v2(i1) /= l2
Next
' 収束判定
' 収束した場合
If Math.Abs((l2-l1)/l1) < eps
Dim k1 As Integer = -1
Dim ii As Integer = 0
Do While ii < n and k1 < 0
If Math.Abs(v2(ii)) > 0.001
k1 = ii
If v2(k1)*v1(k1) < 0.0
l2 = -l2
End If
End If
ii += 1
Loop
k = ct
r(i) = l2
For i1 As Integer = 0 To n-1
v(i,i1) = v2(i1)
Next
If i = n-1
ind = i + 1
Else
For i1 As Integer = 0 To n-1
For i2 As Integer = 0 To n-1
C(i1,i2) = 0.0
For i3 As Integer = 0 To n-1
Dim x As Double = A(i1,i3) - l2
If i1 <> i3
x = A(i1,i3)
End If
C(i1,i2) += x * B(i3,i2)
Next
Next
Next
For i1 As Integer = 0 To n-1
For i2 As Integer = 0 To n-1
B(i1,i2) = C(i1,i2)
Next
Next
For i1 As Integer = 0 To n-1
v1(i1) = 0.0
For i2 As Integer = 0 To n-1
v1(i1) += B(i1,i2) * v0(i2)
Next
Next
For i1 As Integer = 0 To n-1
v0(i1) = v1(i1)
Next
ind = power(i+1, n, m, ct, eps, A, B, C, r, v, v0, v1, v2)
End If
' 収束しない場合
Else
For i1 As Integer = 0 To n-1
v1(i1) = v2(i1)
Next
l1 = l2
k += 1
End If
Loop
Return ind
End Function
''''''''''''''''''''''''''''''''''''''''''
' 線形連立方程式を解く（逆行列を求める） '
' w : 方程式の左辺及び右辺 '
' n : 方程式の数 '
' m : 方程式の右辺の列の数 '
' eps : 逆行列の存在を判定する規準 '
' return : =0 : 正常 '
' =1 : 逆行列が存在しない '
''''''''''''''''''''''''''''''''''''''''''
Function gauss(w(,) As Double, n As Integer, m As Integer, eps As Double) As Integer
Dim ind As Integer = 0
Dim nm As Integer = n + m
Dim i1 As Integer = 0
Do While i1 < n and ind = 0
Dim y1 As Double = 0.0
Dim m1 As Integer = i1 + 1
Dim m2 As Integer = 0
' ピボット要素の選択
For i2 As Integer = i1 To n-1
Dim y2 As Double = Math.Abs(w(i2,i1))
If y1 < y2
y1 = y2
m2 = i2
End If
Next
' 逆行列が存在しない
If y1 < eps
ind = 1
' 逆行列が存在する
Else
' 行の入れ替え
For i2 As Integer = i1 To nm-1
y1 = w(i1,i2)
w(i1,i2) = w(m2,i2)
w(m2,i2) = y1
Next
' 掃き出し操作
y1 = 1.0 / w(i1,i1)
For i2 As Integer = m1 To nm-1
w(i1,i2) *= y1
Next
For i2 As Integer = 0 To n-1
If i2 <> i1
For i3 As Integer = m1 To nm-1
w(i2,i3) -= w(i2,i1) * w(i1,i3)
Next
End If
Next
End If
i1 += 1
Loop
Return ind
End Function
End Module
/*
---------データ例（コメント部分を除いて下さい）---------
2 2 100 // 各組の変数の数（r と s， r ≦ s）とデータの数（n）
66 22 44 31 // x1, x2, x3, x4
25 74 17 81
50 23 53 71
25 57 19 81
74 47 64 47
39 33 48 46
14 22 9 69
67 60 49 26
42 40 77 65
11 80 0 86
32 0 43 74
68 69 44 68
24 49 9 71
42 74 28 46
60 58 73 28
36 37 33 68
24 44 19 83
30 40 31 50
55 40 60 49
63 47 94 41
72 30 100 45
19 22 13 75
43 39 43 34
90 83 92 31
51 77 52 82
53 70 34 31
28 51 53 44
40 62 42 79
31 48 22 68
57 29 51 30
64 89 57 42
49 82 72 29
53 31 55 43
79 52 70 10
45 19 43 57
35 34 34 89
4 69 0 100
49 49 66 66
92 82 97 6
5 89 0 100
65 26 83 28
56 36 64 38
48 50 25 22
30 30 15 55
40 65 38 42
14 67 9 67
84 96 90 8
53 64 51 54
50 89 60 52
76 41 68 9
49 40 53 49
78 66 66 17
76 58 90 29
41 15 40 49
63 60 55 33
40 36 49 67
78 54 71 18
62 72 69 12
64 47 42 53
56 64 9 15
77 35 56 25
44 12 46 87
80 9 56 19
36 21 52 78
48 63 64 48
43 61 50 47
58 23 28 50
90 12 100 0
13 33 11 77
67 44 48 28
75 45 68 17
81 22 89 9
46 45 59 55
56 49 64 55
65 62 72 27
34 49 29 77
45 33 60 63
20 45 14 99
33 38 26 87
44 51 69 52
64 57 64 48
44 64 51 28
63 48 56 11
29 39 33 84
40 48 51 54
40 38 26 62
68 46 61 26
58 45 68 48
64 44 77 63
59 62 44 66
81 53 93 19
23 34 12 68
51 35 55 46
74 70 84 17
42 33 56 44
46 31 46 53
33 57 38 63
40 24 20 42
53 36 60 31
0 34 0 100
*/
| 情報学部 | 菅沼ホーム | 目次 | 索引 |