| 情報学部 | 菅沼ホーム | 目次 | 索引 |
/****************************/
/* 主成分分析 */
/* coded by Y.Suganuma */
/****************************/
#include <stdio.h>
#include <math.h>
int principal(int, int, double **, double *, double **, double, int);
int Jacobi(int, int, double, double **, double **, double **, double **, double **);
int main()
{
double **x, *r, **a;
int i1, i2, n, p, sw;
scanf("%d %d", &p, &n); // 変数の数とデータの数
r = new double [p];
x = new double * [p];
a = new double * [p];
for (i1 = 0; i1 < p; i1++) {
x[i1] = new double [n];
a[i1] = new double [p];
}
for (i1 = 0; i1 < n; i1++) { // データ
for (i2 = 0; i2 < p; i2++)
scanf("%lf", &x[i2][i1]);
}
sw = principal(p, n, x, r, a, 1.0e-10, 200);
if (sw == 0) {
for (i1 = 0; i1 < p; i1++) {
printf("主成分 %f", r[i1]);
printf(" 係数");
for (i2 = 0; i2 < p; i2++)
printf(" %f", a[i1][i2]);
printf("\n");
}
}
else
printf("***error 解を求めることができませんでした\n");
for (i1 = 0; i1 < p; i1++) {
delete [] x[i1];
delete [] a[i1];
}
delete [] x;
delete [] a;
delete [] r;
return 0;
}
/***********************************/
/* 主成分分析 */
/* p : 変数の数 */
/* n : データの数 */
/* x : データ */
/* r : 分散（主成分） */
/* a : 係数 */
/* eps : 収束性を判定する規準 */
/* ct : 最大繰り返し回数 */
/* return : =0 : 正常 */
/* =1 : エラー */
/***********************************/
int principal(int p, int n, double **x, double *r, double **a, double eps, int ct)
{
double **A1, **A2, **C, mean, **X1, **X2, s2;
int i1, i2, i3, sw = 0;
// 領域の確保
C = new double * [p];
A1 = new double * [p];
A2 = new double * [p];
X1 = new double * [p];
X2 = new double * [p];
for (i1 = 0; i1 < p; i1++) {
C[i1] = new double [p];
A1[i1] = new double [p];
A2[i1] = new double [p];
X1[i1] = new double [p];
X2[i1] = new double [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;
}
}
// 固有値と固有ベクトルの計算（ヤコビ法）
sw = Jacobi(p, ct, eps, C, A1, A2, X1, X2);
if (sw == 0) {
for (i1 = 0; i1 < p; i1++) {
r[i1] = A1[i1][i1];
for (i2 = 0; i2 < p; i2++)
a[i1][i2] = X1[i2][i1];
}
}
// 領域の解放
for (i1 = 0; i1 < p; i1++) {
delete [] C[i1];
delete [] A1[i1];
delete [] A2[i1];
delete [] X1[i1];
delete [] X2[i1];
}
delete [] C;
delete [] A1;
delete [] A2;
delete [] X1;
delete [] X2;
return sw;
}
/*************************************************************/
/* 実対称行列の固有値・固有ベクトル（ヤコビ法） */
/* n : 次数 */
/* ct : 最大繰り返し回数 */
/* eps : 収束判定条件 */
/* A : 対象とする行列 */
/* A1, A2 : 作業域（nxnの行列），A1の対角要素が固有値 */
/* X1, X2 : 作業域（nxnの行列），X1の各列が固有ベクトル */
/* return : =0 : 正常 */
/* =1 : 収束せず */
/* coded by Y.Suganuma */
/*************************************************************/
#include <math.h>
int Jacobi(int n, int ct, double eps, double **A, double **A1, double **A2,
double **X1, double **X2)
{
double max, s, t, v, sn, cs;
int i1, i2, 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 (fabs(A1[i1][i2]) > max) {
max = fabs(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 = fabs(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;
}
/*
---------データ例（コメント部分を除いて下さい）---------
4 100 // 変数の数（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
*/
/****************************/
/* 主成分分析 */
/* coded by Y.Suganuma */
/****************************/
import java.io.*;
import java.util.StringTokenizer;
public class Test {
public static void main(String args[]) throws IOException
{
double x[][], r[], a[][];
int i1, i2, n, p, sw;
StringTokenizer str;
BufferedReader in = new BufferedReader(new InputStreamReader(System.in));
// 変数の数とデータの数
str = new StringTokenizer(in.readLine(), " ");
p = Integer.parseInt(str.nextToken());
n = Integer.parseInt(str.nextToken());
r = new double [p];
x = new double [p][n];
a = new double [p][p];
// データ
for (i1 = 0; i1 < n; i1++) {
str = new StringTokenizer(in.readLine(), " ");
for (i2 = 0; i2 < p; i2++)
x[i2][i1] = Double.parseDouble(str.nextToken());
}
sw = principal(p, n, x, r, a, 1.0e-10, 200);
if (sw == 0) {
for (i1 = 0; i1 < p; i1++) {
System.out.print("主成分 " + r[i1]);
System.out.print(" 係数");
for (i2 = 0; i2 < p; i2++)
System.out.print(" " + a[i1][i2]);
System.out.println();
}
}
else
System.out.println("***error 解を求めることができませんでした");
}
/***********************************/
/* 主成分分析 */
/* p : 変数の数 */
/* n : データの数 */
/* x : データ */
/* r : 分散（主成分） */
/* a : 係数 */
/* eps : 収束性を判定する規準 */
/* ct : 最大繰り返し回数 */
/* return : =0 : 正常 */
/* =1 : エラー */
/***********************************/
static int principal(int p, int n, double x[][], double r[], double a[][], double eps, int ct)
{
double A1[][], A2[][], C[][], mean, X1[][], X2[][], s2;
int i1, i2, i3, sw = 0;
// 領域の確保
C = new double [p][p];
A1 = new double [p][p];
A2 = new double [p][p];
X1 = new double [p][p];
X2 = new double [p][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 = Math.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;
}
}
// 固有値と固有ベクトルの計算（ヤコビ法）
sw = Jacobi(p, ct, eps, C, A1, A2, X1, X2);
if (sw == 0) {
for (i1 = 0; i1 < p; i1++) {
r[i1] = A1[i1][i1];
for (i2 = 0; i2 < p; i2++)
a[i1][i2] = X1[i2][i1];
}
}
return sw;
}
/*************************************************************/
/* 実対称行列の固有値・固有ベクトル（ヤコビ法） */
/* n : 次数 */
/* ct : 最大繰り返し回数 */
/* eps : 収束判定条件 */
/* A : 対象とする行列 */
/* A1, A2 : 作業域（nxnの行列），A1の対角要素が固有値 */
/* X1, X2 : 作業域（nxnの行列），X1の各列が固有ベクトル */
/* return : =0 : 正常 */
/* =1 : 収束せず */
/* coded by Y.Suganuma */
/*************************************************************/
static int Jacobi(int n, int ct, double eps, double A[][], double A1[][], double A2[][],
double X1[][], double X2[][])
{
double max, s, t, v, sn, cs;
int i1, i2, 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 (Math.abs(A1[i1][i2]) > max) {
max = Math.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 = Math.abs(t) / Math.sqrt(s * s + t * t);
sn = Math.sqrt(0.5 * (1.0 - v));
if (s*t < 0.0)
sn = -sn;
cs = Math.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;
}
}
/*
---------データ例（コメント部分を除いて下さい）---------
4 100 // 変数の数（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
*/
<!DOCTYPE HTML>
<HTML>
<HEAD>
<TITLE>主成分分析</TITLE>
<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=utf-8">
<SCRIPT TYPE="text/javascript">
function main()
{
// データ
let p = parseInt(document.getElementById("p").value);
let n = parseInt(document.getElementById("n").value);
let r = new Array();
let x = new Array();
let a = new Array();
for (let i1 = 0; i1 < p; i1++) {
x[i1] = new Array();
a[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 < p; i2++) {
x[i2][i1] = parseFloat(ss[k]);
k++;
}
}
// 計算
let sw = principal(p, n, x, r, a, 1.0e-10, 200);
if (sw > 0)
document.getElementById("ans").value = "収束しません";
else {
let str = "";
for (let i1 = 0; i1 < p; i1++) {
str = str + "主成分 " + r[i1] + "\n";
str = str + " 係数";
for (let i2 = 0; i2 < p; i2++)
str = str + " " + a[i1][i2];
str = str + "\n";
}
document.getElementById("ans").value = str;
}
}
/***********************************/
/* 主成分分析 */
/* p : 変数の数 */
/* n : データの数 */
/* x : データ */
/* r : 分散（主成分） */
/* a : 係数 */
/* eps : 収束性を判定する規準 */
/* ct : 最大繰り返し回数 */
/* return : =0 : 正常 */
/* =1 : エラー */
/***********************************/
function principal(p, n, x, r, a, eps, ct)
{
let mean;
let s2;
let i1;
let i2;
let i3;
let sw = 0;
// 領域の確保
let C = new Array();
let A1 = new Array();
let A2 = new Array();
let X1 = new Array();
let X2 = new Array();
for (let i1 = 0; i1 < p; i1++) {
C[i1] = new Array();
A1[i1] = new Array();
A2[i1] = new Array();
X1[i1] = new Array();
X2[i1] = new Array();
}
// データの基準化
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 = Math.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;
}
}
// 固有値と固有ベクトルの計算（ヤコビ法）
sw = Jacobi(p, ct, eps, C, A1, A2, X1, X2);
if (sw == 0) {
for (i1 = 0; i1 < p; i1++) {
r[i1] = A1[i1][i1];
for (i2 = 0; i2 < p; i2++)
a[i1][i2] = X1[i2][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)
{
let max;
let s;
let t;
let v;
let sn;
let cs;
let i1;
let i2;
let k = 0;
let ind = 1;
let p = 0;
let 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 (Math.abs(A1[i1][i2]) > max) {
max = Math.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 = Math.abs(t) / Math.sqrt(s * s + t * t);
sn = Math.sqrt(0.5 * (1.0 - v));
if (s*t < 0.0)
sn = -sn;
cs = Math.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;
}
</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">
変数の数：<INPUT ID="p" STYLE="font-size: 100%" TYPE="text" SIZE="2" VALUE="4">
データの数：<INPUT ID="n" STYLE="font-size: 100%" TYPE="text" SIZE="3" 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", $p, $n); // 変数の数とデータの数
$r = 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 = principal($p, $n, $x, $r, $a, 1.0e-10, 200);
if ($sw == 0) {
for ($i1 = 0; $i1 < $p; $i1++) {
printf("主成分 %f", $r[$i1]);
printf(" 係数");
for ($i2 = 0; $i2 < $p; $i2++)
printf(" %f", $a[$i1][$i2]);
printf("\n");
}
}
else
printf("***error 解を求めることができませんでした\n");
/***********************************/
/* 主成分分析 */
/* p : 変数の数 */
/* n : データの数 */
/* x : データ */
/* r : 分散（主成分） */
/* a : 係数 */
/* eps : 収束性を判定する規準 */
/* ct : 最大繰り返し回数 */
/* return : =0 : 正常 */
/* =1 : エラー */
/***********************************/
function principal($p, $n, $x, &$r, &$a, $eps, $ct)
{
$sw = 0;
// 領域の確保
$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;
}
}
// 固有値と固有ベクトルの計算（ヤコビ法）
$sw = Jacobi($p, $ct, $eps, $C, $A1, $A2, $X1, $X2);
if ($sw == 0) {
for ($i1 = 0; $i1 < $p; $i1++) {
$r[$i1] = $A1[$i1][$i1];
for ($i2 = 0; $i2 < $p; $i2++)
$a[$i1][$i2] = $X1[$i2][$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;
}
/*
---------データ例（コメント部分を除いて下さい）---------
4 100 // 変数の数（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
*/
?>
############################
# 主成分分析
# coded by Y.Suganuma
############################
#************************************************************/
# 実対称行列の固有値・固有ベクトル（ヤコビ法） */
# n : 次数 */
# ct : 最大繰り返し回数 */
# eps : 収束判定条件 */
# a : 対象とする行列 */
# a1, a2 : 作業域（nxnの行列），a1の対角要素が固有値 */
# x1, x2 : 作業域（nxnの行列），x1の各列が固有ベクトル */
# return : =0 : 正常 */
# =1 : 収束せず */
# coded by Y.Suganuma */
#************************************************************/
def Jacobi(n, ct, eps, a, a1, a2, x1, x2)
# 初期設定
k = 0
ind = 1
p = 0
q = 0
for i1 in 0 ... n
for i2 in 0 ... n
a1[i1][i2] = a[i1][i2]
x1[i1][i2] = 0.0
end
x1[i1][i1] = 1.0
end
# 計算
while ind > 0 && k < ct
# 最大要素の探索
max = 0.0
for i1 in 0 ... n
for i2 in 0 ... n
if i2 != i1
if a1[i1][i2].abs() > max
max = a1[i1][i2].abs()
p = i1
q = i2
end
end
end
end
# 収束判定
# 収束した
if max < eps
ind = 0
# 収束しない
else
# 準備
s = -a1[p][q]
t = 0.5 * (a1[p][p] - a1[q][q])
v = t.abs() / Math.sqrt(s * s + t * t)
sn = Math.sqrt(0.5 * (1.0 - v))
if s*t < 0.0
sn = -sn
end
cs = Math.sqrt(1.0 - sn * sn)
# akの計算
for i1 in 0 ... n
if i1 == p
for i2 in 0 ... n
if i2 == p
a2[p][p] = a1[p][p] * cs * cs + a1[q][q] * sn * sn -
2.0 * a1[p][q] * sn * cs
elsif i2 == q
a2[p][q] = 0.0
else
a2[p][i2] = a1[p][i2] * cs - a1[q][i2] * sn
end
end
elsif i1 == q
for i2 in 0 ... n
if (i2 == q)
a2[q][q] = a1[p][p] * sn * sn + a1[q][q] * cs * cs +
2.0 * a1[p][q] * sn * cs
elsif i2 == p
a2[q][p] = 0.0
else
a2[q][i2] = a1[q][i2] * cs + a1[p][i2] * sn
end
end
else
for i2 in 0 ... n
if i2 == p
a2[i1][p] = a1[i1][p] * cs - a1[i1][q] * sn
elsif i2 == q
a2[i1][q] = a1[i1][q] * cs + a1[i1][p] * sn
else
a2[i1][i2] = a1[i1][i2]
end
end
end
end
# xkの計算
for i1 in 0 ... n
for i2 in 0 ... n
if i2 == p
x2[i1][p] = x1[i1][p] * cs - x1[i1][q] * sn
elsif i2 == q
x2[i1][q] = x1[i1][q] * cs + x1[i1][p] * sn
else
x2[i1][i2] = x1[i1][i2]
end
end
end
# 次のステップへ
k += 1
for i1 in 0 ... n
for i2 in 0 ... n
a1[i1][i2] = a2[i1][i2]
x1[i1][i2] = x2[i1][i2]
end
end
end
end
return ind
end
###################################
# 主成分分析
# p : 変数の数
# n : データの数
# x : データ
# r : 分散（主成分）
# a : 係数
# eps : 収束性を判定する規準
# ct : 最大繰り返し回数
# return : =0 : 正常
# =1 : エラー
# coded by Y.Suganuma
###################################
def principal(p, n, x, r, a, eps, ct)
sw = 0
# 領域の確保
c = Array.new(p)
a1 = Array.new(p)
a2 = Array.new(p)
x1 = Array.new(p)
x2 = Array.new(p)
for i1 in 0 ... p
c[i1] = Array.new(p)
a1[i1] = Array.new(p)
a2[i1] = Array.new(p)
x1[i1] = Array.new(p)
x2[i1] = Array.new(p)
end
# データの基準化
for i1 in 0 ... p
mean = 0.0
s2 = 0.0
for i2 in 0 ... n
mean += x[i1][i2]
s2 += x[i1][i2] * x[i1][i2]
end
mean /= n
s2 /= n
s2 = n * (s2 - mean * mean) / (n - 1)
s2 = Math.sqrt(s2)
for i2 in 0 ... n
x[i1][i2] = (x[i1][i2] - mean) / s2
end
end
# 分散強分散行列の計算
for i1 in 0 ... p
for i2 in 0 ... p
s2 = 0.0
for i3 in 0 ... n
s2 += x[i1][i3] * x[i2][i3]
end
s2 /= (n - 1)
c[i1][i2] = s2
if i1 != i2
c[i2][i1] = s2
end
end
end
# 固有値と固有ベクトルの計算（ヤコビ法）
sw = Jacobi(p, ct, eps, c, a1, a2, x1, x2)
if sw == 0
for i1 in 0 ... p
r[i1] = a1[i1][i1]
for i2 in 0 ... p
a[i1][i2] = x1[i2][i1]
end
end
end
return sw
end
ss = gets().split(" ")
p = Integer(ss[0]) # 変数の数
n = Integer(ss[1]) # データの数
r = Array.new(p)
x = Array.new(p)
a = Array.new(p)
for i1 in 0 ... p
x[i1] = Array.new(n)
a[i1] = Array.new(p)
end
for i1 in 0 ... n # データ
ss = gets().split(" ")
for i2 in 0 ... p
x[i2][i1] = Float(ss[i2])
end
end
sw = principal(p, n, x, r, a, 1.0e-10, 200)
if sw == 0
for i1 in 0 ... p
print("主成分 " + String(r[i1]))
print(" 係数")
for i2 in 0 ... p
print(" " + String(a[i1][i2]))
end
print("\n")
end
else
print("***error 解を求めることができませんでした")
end
=begin
---------データ例（コメント部分を除いて下さい）---------
4 100 # 変数の数（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
=end
# -*- coding: UTF-8 -*-
import numpy as np
import sys
from math import *
############################################
# 実対称行列の固有値・固有ベクトル（ヤコビ法）
# n : 次数
# ct : 最大繰り返し回数
# eps : 収束判定条件
# A : 対象とする行列
# A1, A2 : 作業域（nxnの行列），A1の対角要素が固有値
# X1, X2 : 作業域（nxnの行列），X1の各列が固有ベクトル
# return : =0 : 正常
# =1 : 収束せず
# coded by Y.Suganuma
############################################
def Jacobi(n, ct, eps, A, A1, A2, X1, X2) :
# 初期設定
k = 0
ind = 1
p = 0
q = 0
for i1 in range(0, n) :
for i2 in range(0, n) :
A1[i1][i2] = A[i1][i2]
X1[i1][i2] = 0.0
X1[i1][i1] = 1.0
# 計算
while ind > 0 and k < ct :
# 最大要素の探索
max = 0.0
for i1 in range(0, n) :
for i2 in range(0, n) :
if i2 != i1 :
if abs(A1[i1][i2]) > max :
max = fabs(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 = fabs(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 in range(0, n) :
if i1 == p :
for i2 in range(0, n) :
if i2 == p :
A2[p][p] = A1[p][p] * cs * cs + A1[q][q] * sn * sn - 2.0 * A1[p][q] * sn * cs
elif i2 == q :
A2[p][q] = 0.0
else :
A2[p][i2] = A1[p][i2] * cs - A1[q][i2] * sn
elif i1 == q :
for i2 in range(0, n) :
if i2 == q :
A2[q][q] = A1[p][p] * sn * sn + A1[q][q] * cs * cs + 2.0 * A1[p][q] * sn * cs
elif i2 == p :
A2[q][p] = 0.0
else :
A2[q][i2] = A1[q][i2] * cs + A1[p][i2] * sn
else :
for i2 in range(0, n) :
if i2 == p :
A2[i1][p] = A1[i1][p] * cs - A1[i1][q] * sn
elif i2 == q :
A2[i1][q] = A1[i1][q] * cs + A1[i1][p] * sn
else :
A2[i1][i2] = A1[i1][i2]
# Xkの計算
for i1 in range(0, n) :
for i2 in range(0, n) :
if i2 == p :
X2[i1][p] = X1[i1][p] * cs - X1[i1][q] * sn
elif i2 == q :
X2[i1][q] = X1[i1][q] * cs + X1[i1][p] * sn
else :
X2[i1][i2] = X1[i1][i2]
# 次のステップへ
k += 1
for i1 in range(0, n) :
for i2 in range(0, n) :
A1[i1][i2] = A2[i1][i2]
X1[i1][i2] = X2[i1][i2]
return ind
###################################
# 主成分分析
# p : 変数の数
# n : データの数
# x : データ
# r : 分散（主成分）
# a : 係数
# eps : 収束性を判定する規準
# ct : 最大繰り返し回数
# return : =0 : 正常
# =1 : エラー
# coded by Y.Suganuma
###################################
def principal(p, n, x, r, a, eps, ct) :
sw = 0
# 領域の確保
C = np.empty((p, p), np.float)
A1 = np.empty((p, p), np.float)
A2 = np.empty((p, p), np.float)
X1 = np.empty((p, p), np.float)
X2 = np.empty((p, p), np.float)
# データの基準化
for i1 in range(0, p) :
mean = 0.0
s2 = 0.0
for i2 in range(0, n) :
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 in range(0, n) :
x[i1][i2] = (x[i1][i2] - mean) / s2
# 分散強分散行列の計算
for i1 in range(0, p) :
for i2 in range(0, p) :
s2 = 0.0
for i3 in range(0, n) :
s2 += x[i1][i3] * x[i2][i3]
s2 /= (n - 1)
C[i1][i2] = s2
if i1 != i2 :
C[i2][i1] = s2
# 固有値と固有ベクトルの計算（ヤコビ法）
sw = Jacobi(p, ct, eps, C, A1, A2, X1, X2)
if sw == 0 :
for i1 in range(0, p) :
r[i1] = A1[i1][i1]
for i2 in range(0, p) :
a[i1][i2] = X1[i2][i1]
return sw
############################
# 主成分分析
# coded by Y.Suganuma
############################
line = sys.stdin.readline()
ss = line.split()
p = int(ss[0]) # 変数の数
n = int(ss[1]) # データの数
r = np.empty(p, np.float)
x = np.empty((p, n), np.float)
a = np.empty((p, p), np.float)
for i1 in range(0, n) : # データ
line = sys.stdin.readline()
ss = line.split()
for i2 in range(0, p) :
x[i2][i1] = float(ss[i2])
sw = principal(p, n, x, r, a, 1.0e-10, 200)
if sw == 0 :
for i1 in range(0, p) :
print("主成分 " + str(r[i1]), end="")
print(" 係数", end="")
for i2 in range(0, p) :
print(" " + str(a[i1][i2]), end="")
print()
else :
print("***error 解を求めることができませんでした")
"""
---------データ例（コメント部分を除いて下さい）---------
4 100 # 変数の数（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
"""
/****************************/
/* 主成分分析 */
/* 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 p = int.Parse(str[0]);
int n = int.Parse(str[1]);
double[] r = new double [p];
double[][] x = new double [p][];
double[][] a = new double [p][];
for (int i1 = 0; i1 < p; i1++) {
x[i1] = new double [n];
a[i1] = new double [p];
}
// データ
for (int i1 = 0; i1 < n; i1++) {
str = Console.ReadLine().Split(charSep, StringSplitOptions.RemoveEmptyEntries);
for (int i2 = 0; i2 < p; i2++)
x[i2][i1] = double.Parse(str[i2]);
}
int sw = principal(p, n, x, r, a, 1.0e-10, 200);
if (sw == 0) {
for (int i1 = 0; i1 < p; i1++) {
Console.Write("主成分 " + r[i1]);
Console.Write(" 係数");
for (int i2 = 0; i2 < p; i2++)
Console.Write(" " + a[i1][i2]);
Console.WriteLine();
}
}
else
Console.WriteLine("***error 解を求めることができませんでした");
}
/***********************************/
/* 主成分分析 */
/* p : 変数の数 */
/* n : データの数 */
/* x : データ */
/* r : 分散（主成分） */
/* a : 係数 */
/* eps : 収束性を判定する規準 */
/* ct : 最大繰り返し回数 */
/* return : =0 : 正常 */
/* =1 : エラー */
/***********************************/
int principal(int p, int n, double[][] x, double[] r, double[][] a, double eps, int ct)
{
// 領域の確保
double[][] C = new double [p][];
double[][] A1 = new double [p][];
double[][] A2 = new double [p][];
double[][] X1 = new double [p][];
double[][] X2 = new double [p][];
for (int i1 = 0; i1 < p; i1++) {
C[i1] = new double [p];
A1[i1] = new double [p];
A2[i1] = new double [p];
X1[i1] = new double [p];
X2[i1] = new double [p];
}
// データの基準化
for (int i1 = 0; i1 < p; i1++) {
double mean = 0.0;
double s2 = 0.0;
for (int 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 = Math.Sqrt(s2);
for (int i2 = 0; i2 < n; i2++)
x[i1][i2] = (x[i1][i2] - mean) / s2;
}
// 分散強分散行列の計算
for (int i1 = 0; i1 < p; i1++) {
for (int i2 = i1; i2 < p; i2++) {
double s2 = 0.0;
for (int 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;
}
}
// 固有値と固有ベクトルの計算（ヤコビ法）
int sw = Jacobi(p, ct, eps, C, A1, A2, X1, X2);
if (sw == 0) {
for (int i1 = 0; i1 < p; i1++) {
r[i1] = A1[i1][i1];
for (int i2 = 0; i2 < p; i2++)
a[i1][i2] = X1[i2][i1];
}
}
return sw;
}
/*************************************************************/
/* 実対称行列の固有値・固有ベクトル（ヤコビ法） */
/* n : 次数 */
/* ct : 最大繰り返し回数 */
/* eps : 収束判定条件 */
/* A : 対象とする行列 */
/* A1, A2 : 作業域（nxnの行列），A1の対角要素が固有値 */
/* X1, X2 : 作業域（nxnの行列），X1の各列が固有ベクトル */
/* return : =0 : 正常 */
/* =1 : 収束せず */
/* coded by Y.Suganuma */
/*************************************************************/
int Jacobi(int n, int ct, double eps, double[][] A, double[][] A1, double[][] A2,
double[][] X1, double[][] X2)
{
// 初期設定
for (int i1 = 0; i1 < n; i1++) {
for (int i2 = 0; i2 < n; i2++) {
A1[i1][i2] = A[i1][i2];
X1[i1][i2] = 0.0;
}
X1[i1][i1] = 1.0;
}
// 計算
int k = 0, ind = 1;
while (ind > 0 && k < ct) {
// 最大要素の探索
double max = 0.0;
int p = 0, q = 0;
for (int i1 = 0; i1 < n; i1++) {
for (int i2 = 0; i2 < n; i2++) {
if (i2 != i1) {
if (Math.Abs(A1[i1][i2]) > max) {
max = Math.Abs(A1[i1][i2]);
p = i1;
q = i2;
}
}
}
}
// 収束判定
// 収束した
if (max < eps)
ind = 0;
// 収束しない
else {
// 準備
double s = -A1[p][q];
double t = 0.5 * (A1[p][p] - A1[q][q]);
double v = Math.Abs(t) / Math.Sqrt(s * s + t * t);
double sn = Math.Sqrt(0.5 * (1.0 - v));
if (s*t < 0.0)
sn = -sn;
double cs = Math.Sqrt(1.0 - sn * sn);
// Akの計算
for (int i1 = 0; i1 < n; i1++) {
if (i1 == p) {
for (int 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 (int 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 (int 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 (int i1 = 0; i1 < n; i1++) {
for (int 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 (int i1 = 0; i1 < n; i1++) {
for (int i2 = 0; i2 < n; i2++) {
A1[i1][i2] = A2[i1][i2];
X1[i1][i2] = X2[i1][i2];
}
}
}
}
return ind;
}
}
/*
---------データ例（コメント部分を除いて下さい）---------
4 100 // 変数の数（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
*/
'**************************'
' 主成分分析 '
' 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 p As Integer = Integer.Parse(str(0))
Dim n As Integer = Integer.Parse(str(1))
Dim r(p) As Double
Dim x(p,n) As Double
Dim a(p,p) As Double
' データ
For i1 As Integer = 0 To n-1
str = MS.Split(Console.ReadLine().Trim())
For i2 As Integer = 0 To p-1
x(i2,i1) = Double.Parse(str(i2))
Next
Next
Dim sw As Integer = principal(p, n, x, r, a, 1.0e-10, 200)
If sw = 0
For i1 As Integer = 0 To p-1
Console.Write("主成分 " & r(i1))
Console.Write(" 係数")
For i2 As Integer = 0 To p-1
Console.Write(" " & a(i1,i2))
Next
Console.WriteLine()
Next
Else
Console.WriteLine("***error 解を求めることができませんでした")
End If
End Sub
'*********************************'
' 主成分分析 '
' p : 変数の数 '
' n : データの数 '
' x : データ '
' r : 分散（主成分） '
' a : 係数 '
' eps : 収束性を判定する規準 '
' ct : 最大繰り返し回数 '
' return : =0 : 正常 '
' =1 : エラー '
'*********************************'
Function principal(p As Integer, n As Integer, x(,) As Double, r() As Double,
a(,) As Double, eps As Double, ct As Integer)
' 領域の確保
Dim C(p,p) As Double
Dim A1(p,p) As Double
Dim A2(p,p) As Double
Dim X1(p,p) As Double
Dim X2(p,p) As Double
' データの基準化
For i1 As Integer = 0 To p-1
Dim mean As Double = 0.0
Dim s2 As Double = 0.0
For i2 As Integer = 0 To n-1
mean += x(i1,i2)
s2 += x(i1,i2) * x(i1,i2)
Next
mean /= n
s2 /= n
s2 = n * (s2 - mean * mean) / (n - 1)
s2 = Math.Sqrt(s2)
For i2 As Integer = 0 To n-1
x(i1,i2) = (x(i1,i2) - mean) / s2
Next
Next
' 分散強分散行列の計算
For i1 As Integer = 0 To p-1
For i2 As Integer = i1 To p-1
Dim s2 As Double = 0.0
For i3 As Integer = 0 To n-1
s2 += x(i1,i3) * x(i2,i3)
Next
s2 /= (n - 1)
C(i1,i2) = s2
If i1 <> i2
C(i2,i1) = s2
End If
Next
Next
' 固有値と固有ベクトルの計算（ヤコビ法）
Dim sw As Integer = Jacobi(p, ct, eps, C, A1, A2, X1, X2)
If sw = 0
For i1 As Integer = 0 To p-1
r(i1) = A1(i1,i1)
For i2 As Integer = 0 To p-1
a(i1,i2) = X1(i2,i1)
Next
Next
End If
Return sw
End Function
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' 実対称行列の固有値・固有ベクトル（ヤコビ法） '
' n : 次数 '
' ct : 最大繰り返し回数 '
' eps : 収束判定条件 '
' A : 対象とする行列 '
' A1, A2 : 作業域（nxnの行列），A1の対角要素が固有値 '
' X1, X2 : 作業域（nxnの行列），X1の各列が固有ベクトル '
' return : =0 : 正常 '
' =1 : 収束せず '
' coded by Y.Suganuma '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function Jacobi(n As Integer, ct As Integer, eps As Double, A(,) As Double,
A1(,) As Double, A2(,) As Double, X1(,) As Double, X2(,) As Double)
' 初期設定
For i1 As Integer = 0 To n-1
For i2 As Integer = 0 To n-1
A1(i1,i2) = A(i1,i2)
X1(i1,i2) = 0.0
Next
X1(i1,i1) = 1.0
Next
' 計算
Dim k As Integer = 0
Dim ind As Integer = 1
Do While ind > 0 and k < ct
' 最大要素の探索
Dim max As Double = 0.0
Dim p As Integer = 0
Dim q As Integer = 0
For i1 As Integer = 0 To n-1
For i2 As Integer = 0 To n-1
If i2 <> i1
If Math.Abs(A1(i1,i2)) > max
max = Math.Abs(A1(i1,i2))
p = i1
q = i2
End If
End If
Next
Next
' 収束判定
' 収束した
If max < eps
ind = 0
' 収束しない
Else
' 準備
Dim s As Double = -A1(p,q)
Dim t As Double = 0.5 * (A1(p,p) - A1(q,q))
Dim v As Double = Math.Abs(t) / Math.Sqrt(s * s + t * t)
Dim sn As Double = Math.Sqrt(0.5 * (1.0 - v))
If s*t < 0.0
sn = -sn
End If
Dim cs As Double = Math.Sqrt(1.0 - sn * sn)
' Akの計算
For i1 As Integer = 0 To n-1
If i1 = p
For i2 As Integer = 0 To n-1
If i2 = p
A2(p,p) = A1(p,p) * cs * cs + A1(q,q) * sn * sn -
2.0 * A1(p,q) * sn * cs
ElseIf i2 = q
A2(p,q) = 0.0
Else
A2(p,i2) = A1(p,i2) * cs - A1(q,i2) * sn
End If
Next
ElseIf i1 = q
For i2 As Integer = 0 To n-1
If i2 = q
A2(q,q) = A1(p,p) * sn * sn + A1(q,q) * cs * cs +
2.0 * A1(p,q) * sn * cs
ElseIf i2 = p
A2(q,p) = 0.0
Else
A2(q,i2) = A1(q,i2) * cs + A1(p,i2) * sn
End If
Next
Else
For i2 As Integer = 0 To n-1
If i2 = p
A2(i1,p) = A1(i1,p) * cs - A1(i1,q) * sn
ElseIf i2 = q
A2(i1,q) = A1(i1,q) * cs + A1(i1,p) * sn
Else
A2(i1,i2) = A1(i1,i2)
End If
Next
End If
Next
' Xkの計算
For i1 As Integer = 0 To n-1
For i2 As Integer = 0 To n-1
If i2 = p
X2(i1,p) = X1(i1,p) * cs - X1(i1,q) * sn
ElseIf i2 = q
X2(i1,q) = X1(i1,q) * cs + X1(i1,p) * sn
Else
X2(i1,i2) = X1(i1,i2)
End If
Next
Next
' 次のステップへ
k += 1
For i1 As Integer = 0 To n-1
For i2 As Integer = 0 To n-1
A1(i1,i2) = A2(i1,i2)
X1(i1,i2) = X2(i1,i2)
Next
Next
End If
Loop
Return ind
End Function
End Module
/*
---------データ例（コメント部分を除いて下さい）---------
4 100 // 変数の数（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
*/
| 情報学部 | 菅沼ホーム | 目次 | 索引 |