因子分析

<?php

/****************************/
/* 因子分析                 */
/*      coded by Y.Suganuma */
/****************************/

	fscanf(STDIN, "%d %d %d", $m, $p, $n);   // 変数の数とデータの数

	$iv = array($p);
	$x  = array($p);
	$a  = array($p);
	for ($i1 = 0; $i1 < $p; $i1++) {
		$x[$i1] = array($n);
		$a[$i1] = array($p);
	}

	for ($i1 = 0; $i1 < $n; $i1++) {   // データ
		$str       = trim(fgets(STDIN));
		$x[0][$i1] = floatval(strtok($str, " "));
		for ($i2 = 1; $i2 < $p; $i2++)
			$x[$i2][$i1] = floatval(strtok(" "));
	}

	$sw = factor($p, $n, $m, $x, $iv, $a, 1.0e-10, 200);

	if ($sw == 0) {
		for ($i1 = 0; $i1 < $m; $i1++) {
			printf("固有値 %f", $iv[$i1]);
			printf(" 共通因子負荷量");
			for ($i2 = 0; $i2 < $p; $i2++)
				printf(" %f", $a[$i1][$i2]);
			printf("\n");
		}
	}
	else
		printf("***error  解を求めることができませんでした\n");

/***********************************/
/* 因子分析                        */
/*      p : 変数の数               */
/*      n : データの数             */
/*      m : 共通因子負荷量の数     */
/*      x : データ                 */
/*      iv : 固有値                */
/*      a : 係数                   */
/*      eps : 収束性を判定する規準 */
/*      ct : 最大繰り返し回数      */
/*      return : =0 : 正常         */
/*               =1 : エラー       */
/***********************************/
function factor($p, $n, $m, $x, &$iv, &$a, $eps, $ct)
{
	$sw    = 1;
	$count = 0;
					// 領域の確保
	$s  = array($p);
	$h1 = array($p);
	$h2 = array($p);
	$C  = array($p);
	$A1 = array($p);
	$A2 = array($p);
	$X1 = array($p);
	$X2 = array($p);
	for ($i1 = 0; $i1 < $p; $i1++) {
		$C[$i1]  = array($p);
		$A1[$i1] = array($p);
		$A2[$i1] = array($p);
		$X1[$i1] = array($p);
		$X2[$i1] = array($p);
	}
					// データの基準化
	for ($i1 = 0; $i1 < $p; $i1++) {
		$mean = 0.0;
		$s2   = 0.0;
		for ($i2 = 0; $i2 < $n; $i2++) {
			$mean += $x[$i1][$i2];
			$s2   += $x[$i1][$i2] * $x[$i1][$i2];
		}
		$mean /= $n;
		$s2   /= $n;
		$s2    = $n * ($s2 - $mean * $mean) / ($n - 1);
		$s2    = sqrt($s2);
		for ($i2 = 0; $i2 < $n; $i2++)
			$x[$i1][$i2] = ($x[$i1][$i2] - $mean) / $s2;
	}
					// 分散強分散行列の計算
	for ($i1 = 0; $i1 < $p; $i1++) {
		for ($i2 = $i1; $i2 < $p; $i2++) {
			$s2 = 0.0;
			for ($i3 = 0; $i3 < $n; $i3++)
				$s2 += $x[$i1][$i3] * $x[$i2][$i3];
			$s2          /= ($n - 1);
			$C[$i1][$i2]  = $s2;
			if ($i1 != $i2)
				$C[$i2][$i1] = $s2;
		}
	}
					// 共通性の初期値
	for ($i1 = 0; $i1 < $p; $i1++)
		$h1[$i1] = 1.0;
					// 共通因子負荷量の計算
	while ($count < $ct && $sw > 0) {
						// 固有値と固有ベクトルの計算（ヤコビ法）
		$sw = Jacobi($p, $ct, $eps, $C, $A1, $A2, $X1, $X2);

		if ($sw == 0) {
						// 共通性及び共通因子負荷量の計算
			for ($i1 = 0; $i1 < $p; $i1++)
				$s[$i1] = 0;

			for ($i1 = 0; $i1 < $m; $i1++) {
				$max_i = -1;
				$max   = 0.0;
				for ($i2 = 0; $i2 < $p; $i2++) {
					if ($s[$i2] == 0 && ($max_i < 0 || $A1[$i2][$i2] > $max)) {
						$max_i = $i2;
						$max   = $A1[$i2][$i2];
					}
				}
				$s[$max_i]  = 1;
				$iv[$i1]    = $A1[$max_i][$max_i];
				$s2         = sqrt($iv[$i1]);
				for ($i2 = 0; $i2 < $p; $i2++)
					$a[$i1][$i2] = $s2 * $X1[$i2][$max_i];
			}

			for ($i1 = 0; $i1 < $p; $i1++) {
				$h2[$i1] = 0.0;
				for ($i2 = 0; $i2 < $m; $i2++)
					$h2[$i1] += $a[$i2][$i1] * $a[$i2][$i1];
			}

			for ($i1 = 0; $i1 < $m && $sw == 0; $i1++) {
				if (abs($h2[$i1]-$h1[$i1]) > $eps)
					$sw = 1;
			}
						// 収束しない場合
			if ($sw > 0) {
				$count++;
				for ($i1 = 0; $i1 < $p; $i1++) {
					$h1[$i1]     = $h2[$i1];
					$C[$i1][$i1] = $h1[$i1];
				}
			}
		}
	}

	return $sw;
}

/*************************************************************/
/* 実対称行列の固有値・固有ベクトル（ヤコビ法）              */
/*      n : 次数                                             */
/*      ct : 最大繰り返し回数                                */
/*      eps : 収束判定条件                                   */
/*      A : 対象とする行列                                   */
/*      A1, A2 : 作業域（nxnの行列），A1の対角要素が固有値   */
/*      X1, X2 : 作業域（nxnの行列），X1の各列が固有ベクトル */
/*      return : =0 : 正常                                   */
/*               =1 : 収束せず                               */
/*      coded by Y.Suganuma                                  */
/*************************************************************/
function Jacobi($n, $ct, $eps, $A, &$A1, $A2, &$X1, $X2)
{
					// 初期設定
	$k   = 0;
	$ind = 1;
	$p   = 0;
	$q   = 0;

	for ($i1 = 0; $i1 < $n; $i1++) {
		for ($i2 = 0; $i2 < $n; $i2++) {
			$A1[$i1][$i2] = $A[$i1][$i2];
			$X1[$i1][$i2] = 0.0;
		}
		$X1[$i1][$i1] = 1.0; 
	}
					// 計算
	while ($ind > 0 && $k < $ct) {
						// 最大要素の探索
		$max = 0.0;
		for ($i1 = 0; $i1 < $n; $i1++) {
			for ($i2 = 0; $i2 < $n; $i2++) {
				if ($i2 != $i1) {
					if (abs($A1[$i1][$i2]) > $max) {
						$max = abs($A1[$i1][$i2]);
						$p   = $i1;
						$q   = $i2;
					}
				}
			}
		}
						// 収束判定
							// 収束した
		if ($max < $eps)
			$ind = 0;
							// 収束しない
		else {
								// 準備
			$s  = -$A1[$p][$q];
			$t  = 0.5 * ($A1[$p][$p] - $A1[$q][$q]);
			$v  = abs($t) / sqrt($s * $s + $t * $t);
			$sn = sqrt(0.5 * (1.0 - $v));
			if ($s*$t < 0.0)
				$sn = -$sn;
			$cs = sqrt(1.0 - $sn * $sn);
								// Akの計算
			for ($i1 = 0; $i1 < $n; $i1++) {
				if ($i1 == $p) {
					for ($i2 = 0; $i2 < $n; $i2++) {
						if ($i2 == $p)
							$A2[$p][$p] = $A1[$p][$p] * $cs * $cs + $A1[$q][$q] * $sn * $sn -
                                          2.0 * $A1[$p][$q] * $sn * $cs;
						else if ($i2 == $q)
							$A2[$p][$q] = 0.0;
						else
							$A2[$p][$i2] = $A1[$p][$i2] * $cs - $A1[$q][$i2] * $sn;
					}
				}
				else if ($i1 == $q) {
					for ($i2 = 0; $i2 < $n; $i2++) {
						if ($i2 == $q)
							$A2[$q][$q] = $A1[$p][$p] * $sn * $sn + $A1[$q][$q] * $cs * $cs +
                                          2.0 * $A1[$p][$q] * $sn * $cs;
						else if ($i2 == $p)
							$A2[$q][$p] = 0.0;
						else
							$A2[$q][$i2] = $A1[$q][$i2] * $cs + $A1[$p][$i2] * $sn;
					}
				}
				else {
					for ($i2 = 0; $i2 < $n; $i2++) {
						if ($i2 == $p)
							$A2[$i1][$p] = $A1[$i1][$p] * $cs - $A1[$i1][$q] * $sn;
						else if ($i2 == $q)
							$A2[$i1][$q] = $A1[$i1][$q] * $cs + $A1[$i1][$p] * $sn;
						else
							$A2[$i1][$i2] = $A1[$i1][$i2];
					}
				}
			}
								// Xkの計算
			for ($i1 = 0; $i1 < $n; $i1++) {
				for ($i2 = 0; $i2 < $n; $i2++) {
					if ($i2 == $p)
						$X2[$i1][$p] = $X1[$i1][$p] * $cs - $X1[$i1][$q] * $sn;
					else if ($i2 == $q)
						$X2[$i1][$q] = $X1[$i1][$q] * $cs + $X1[$i1][$p] * $sn;
					else
						$X2[$i1][$i2] = $X1[$i1][$i2];
				}
			}
								// 次のステップへ
			$k++;
			for ($i1 = 0; $i1 < $n; $i1++) {
				for ($i2 = 0; $i2 < $n; $i2++) {
					$A1[$i1][$i2] = $A2[$i1][$i2];
					$X1[$i1][$i2] = $X2[$i1][$i2];
				}
			}
		}
	}

	return $ind;
}

/*
---------データ例（コメント部分を除いて下さい）---------
2 4 100   // 共通因子負荷量の数（m），変数の数（p），及び，データの数（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
*/

?>