正準相関分析

<?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
*/

?>