注意点：代入する配列のサイズはn行n列になっていること．メインプログラムでは同じサイズの配列を宣言しておかないとエラーになる．
SUBROUTINE jacobi(a,n,d,v)
INTEGER :: n,np,nrot
DOUBLE PRECISION, INTENT(INOUT) :: a(n,n),d(n),v(n,n)
INTEGER, PARAMETER :: NMAX=500
INTEGER :: i,ip,iq,j
DOUBLE PRECISION :: c,g,h,s,sm,t,tau,theta,tresh,b(NMAX),z(NMAX)
do ip=1,n
do iq=1,n
v(ip,iq)=0.D0
enddo
v(ip,ip)=1.D0
enddo
do ip=1,n
b(ip)=a(ip,ip)
d(ip)=b(ip)
z(ip)=0.D0
enddo
nrot=0
do i=1,50
sm=0.D0
do ip=1,n-1
do iq=ip+1,n
sm=sm+abs(a(ip,iq))
enddo
enddo
if(sm == 0.D0)return
if(i .lt. 4)then
tresh=0.2D0*sm/n**2
else
tresh=0.D0
endif
do ip=1,n-1
do iq=ip+1,n
g=100.D0*dabs(a(ip,iq))
if((i.gt.4).and.(dabs(d(ip))+g.eq.dabs(d(ip))) &
.and.(dabs(d(iq))+g.eq.dabs(d(iq))))then
a(ip,iq)=0.D0
else if(dabs(a(ip,iq)).gt.tresh)then
h=d(iq)-d(ip)
if(dabs(h)+g.eq.dabs(h))then
t=a(ip,iq)/h
else
theta=0.5D0*h/a(ip,iq)
t=1./(dabs(theta)+dsqrt(1.+theta**2))
if(theta.lt.0.) t=-t
endif
c=1.D0/dsqrt(1+t**2)
s=t*c
tau=s/(1.+c)
h=t*a(ip,iq)
z(ip)=z(ip)-h
z(iq)=z(iq)+h
d(ip)=d(ip)-h
d(iq)=d(iq)+h
a(ip,iq)=0.
do j=1,ip-1
g=a(j,ip)
h=a(j,iq)
a(j,ip)=g-s*(h+g*tau)
a(j,iq)=h+s*(g-h*tau)
enddo
do j=ip+1,iq-1
g=a(ip,j)
h=a(j,iq)
a(ip,j)=g-s*(h+g*tau)
a(j,iq)=h+s*(g-h*tau)
enddo
do j=iq+1,n
g=a(ip,j)
h=a(iq,j)
a(ip,j)=g-s*(h+g*tau)
a(iq,j)=h+s*(g-h*tau)
enddo
do j=1,n
g=v(j,ip)
h=v(j,iq)
v(j,ip)=g-s*(h+g*tau)
v(j,iq)=h+s*(g-h*tau)
enddo
nrot=nrot+1
endif
enddo
enddo
do ip=1,n
b(ip)=b(ip)+z(ip)
d(ip)=b(ip)
z(ip)=0.
enddo
enddo
return
END SUBROUTINE jacobi
5.D0 2.D0 2.D0 2.D0
データの種類（国語，数学．．）を \(j\) ，サンプル番号（生徒の番号）を \(i\) で記す．データの種類は \(n\) ，サンプルは \(m\) 個あるものとする．データ種類 \(j\)，サンプル番号 \(i\) の値を\(x_{ij}\)とする．例えば，生徒番号4，数学の成績（練習問題２）は\(x_{43}\)となる．
一般に異なる変量は異なる次元を持つ．そのためまず各変量をおなじように扱えるようにするために標準化という作業をおこなう．これは，平均値がゼロ，分散が1となるような変換をデータに施すことである．変量\(x_{ij}\)は
\begin{equation} z_{ij} = {x_{ij} - \bar{x}_j\over \sigma_{x_j}} \end{equation} と\(z_{ij}\)に変換される．ここで平均値と分散は \begin{eqnarray} \bar{x}_{j} &=& {1\over m}\sum_{i=1}^mx_{ij}\nonumber\\ \sigma_{x_j}^2 &=& {1\over m}\sum_{i=1}^m (x_{ij}-\bar{x}_j)^2 \end{eqnarray} で与えられる． 標準化の結果, 全てのデータ種の平均値が0, 分散が1となる. \begin{eqnarray} \bar{z}_{j} &=& {1\over m}\sum_{i=1}^mz_{ij}=0\nonumber\\ \sigma_{z_j}^2 &=& {1\over m}\sum_{i=1}^m (z_{ij})^2={1\over m}\sum_{i=1}^m\left({x_{ij} - \bar{x}_j\over \sigma_{x_j}}\right)^2=1 \end{eqnarray} 標準化を行った結果，\(i\)を行番号，\(j\)を列番号とする行列\(z_{ij}\)が作られる． \begin{equation} Z= \left( \begin{array}{cccc} z_{11} & z_{12} & \ldots & z_{1n} \\ z_{21} & z_{22} & \ldots & z_{2n} \\ ... & ... & ... & ...\\ z_{m1} & z_{m2} & \ldots & z_{mn} \end{array} \right) \end{equation}
この行列の転置行列を\(Z\)にかけると次の行列Rを得ることができる．
\begin{eqnarray}
R={1\over m}Z^{T}Z&=&{1\over m}
\left(
\begin{array}{cccc}
z_{11} & z_{21} & \ldots & z_{m1} \\
z_{12} & z_{22} & \ldots & z_{m2} \\
... & ... & ... & ...\\
z_{1n} & z_{2n} & \ldots & z_{mn}
\end{array}
\right)
\left(
\begin{array}{cccc}
z_{11} & z_{12} & \ldots & z_{1n} \\
z_{21} & z_{22} & \ldots & z_{2n} \\
... & ... & ... & ...\\
z_{m1} & z_{m2} & \ldots & z_{mn}
\end{array}
\right)\nonumber\\
&=&{1\over m}\left(
\begin{array}{cccc}
\sum_{i=1}^{m} z_{i1}^2 & \sum_{i=1}^{m} z_{i1}z_{i2} & \ldots & \sum_{i=1}^{m} z_{i1}z_{im} \\
\sum_{i=1}^{m} z_{i2}z_{i1} & \sum_{i=1}^{m} z_{i2}^2 & \ldots & \sum_{i=1}^{m} z_{i2}z_{im} \\
... & ... & ... & ...\\
\sum_{i=1}^{m} z_{in}z_{i1} & \sum_{i=1}^{m} z_{in}z_{i2} & \ldots & \sum_{i=1}^{m} z_{in}^2
\end{array}
\right)
\end{eqnarray}
各要素は二つの変量の間の相関係数
\begin{equation}
\rho_{jk}
={1/m\sum_{i=1}^m(x_{ij}-\bar{x}_j)(x_{ik}-\bar{x}_k)\over
\sqrt{\sum_{i=1}^m(x_{ij}-\bar{x}_j)^2/m\sum_{i=1}^m(x_{ik}-\bar{x}_k)^2}/m}
\end{equation}
を与える．この行列を相関行列と呼ぶ．
二つの変量の相関がよいほど相関係数の値の絶対値は1に近 くなる．行列\(R\)の要素はそれぞれ二つの変量の相関係数のため，\(R\)は相関行列と 呼ばれる．
それではでて来た相関係数はどう味わえばよいのか？ それには次の表を用いる. Nはデータの数, rは相関係数である. それぞれの数値は, 無相関の\(N\)個のデータがあった時に\(r\)以上の相関係数が出現する確率（パーセント）である. この数値が小さい程, 二つのデータには相関があると思ってよい. 5以下なら有意, 1以下なら高度に有意であるという.
| r | |||||||||||
| N | 0.0 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.6 | 0.7 | 0.8 | 0.9 | 1.0 |
| 3 | 100 | 94 | 87 | 81 | 74 | 67 | 59 | 51 | 41 | 29 | 0 |
| 6 | 100 | 85 | 70 | 56 | 43 | 31 | 21 | 12 | 6 | 1 | 0 |
| 10 | 100 | 78 | 58 | 40 | 25 | 14 | 7 | 2 | 0.5 | 0 | |
| 20 | 100 | 67 | 40 | 20 | 8 | 2 | 0.5 | 0.1 | 0 | ||
| 50 | 100 | 49 | 16 | 3 | 0.4 | 0 |
無相関な二つの量が \(N\) 組あったとき, 相関係数が\(r\)以上になる確率. 空白は0.05%以下.
以下はテストのデータである （多変量解析入門，永田・棟近共著，サイエンス社 より）．相関行列を生成し，どの科目とどの科目との間に強い相関があるのか確認せよ．
| 生徒番号 i | 国語 x1 | 英語 x2 | 数学 x3 | 理科 x4 |
| 1 | 86 | 79 | 67 | 68 |
| 2 | 71 | 75 | 78 | 84 |
| 3 | 42 | 43 | 39 | 44 |
| 4 | 62 | 58 | 98 | 95 |
| 5 | 96 | 97 | 61 | 63 |
| 6 | 39 | 33 | 45 | 50 |
| 7 | 50 | 53 | 64 | 72 |
| 8 | 78 | 66 | 52 | 47 |
| 9 | 51 | 44 | 76 | 72 |
| 10 | 89 | 92 | 93 | 91 |
宣言文 integer, parameter :: n=4, m=10 double precision :: z(n,m) x(n,m), xsigma(n), xav(n) 成績データの読み込み read(1,*) ( x(j,i), j = 1, n ) 各科目nについて，平均値xav(n)と標準偏差xsigma(n)を求める 式(2) 平均値と標準偏差を用いて，x(n,m)からz(n,m)を作る 式(1), (4) R=ZTZ/mを求め，行列の各要素を出力する 式(5)
\(i\)番目のデータ種に着目したとする．\(z_{ij}\)のヒストグラムを作ったとすると，データは標準化されているので平均がゼロ，分散1の正規分布となる．これはiによらない．つまり，データを標準化すると全てのデータ種について同じ分布となっている．
複数のデータ種の線形結合を作り，この分布がどうなるか考える．全てのデータ種の平均値はゼロなので，線型結合の平均値もゼロ．しかし線型結合の分散は一般に1とは異なってくる．分散が大きい場合と小さい場合とがありえる．分散が小さいということは，その線形結合に対する依存性が小さいということ．逆に分散が大きくなると，その組み合わせによってデータの依存性がよく説明できるということになる．
データの線型結合を
\begin{equation} u_{i}=a_1z_{i1}+a_2z_{i2}+...+a_nz_{in}=\sum_{j=1}^{n}a_jz_{ij} \end{equation}で与える．この線型結合の二乗和（平均値がゼロなので分散と同一）を作る．
\begin{eqnarray} S&=&{1\over m}\sum_{i=1}^m(u_{i})^2\nonumber\\ &=&{1\over m}\sum_{i=1}^m(\sum_{j=1}^{n}a_jz_{ij})^2\nonumber\\ &=&{1\over m}\sum_{i=1}^m(a_1^2z^2_{i1}+a_2^2z^2_{i2}+... +2a_1a_2z_{i1}z_{i2}+a_1a_3z_{i1}z_{i3}+..)\nonumber\\ &=& a_1^2+a_2^2+...+2\rho_{12}a_1a_2+... \end{eqnarray}ここで，各データ種について分散は１，平均はゼロなので，\(z_{ij}\)のjについての二乗和は全て１になることに注意．
また，各iについて平均値はゼロなので，ここで\(\rho_{ik}\)は相関係数となっていることに注意．
\(S\)が「最大」になる\(a_i\)の組み合わせを探す．ただし，\(a_i\)ベクトルの長さは1
\begin{equation} a_1^2+a_2^2+...a_n^2=1 \end{equation}という制約をつけておく必要がある（そうしておかないと\(a\)を大きくすれば\(S\)が大きくなる）．このような条件では，ラグランジュの未定乗数法を用いることで\(a_i\)を見つけることができる．ラグランジュの未定乗数\(\lambda\)を用いて
\begin{equation} f(a_1,a_2,...a_n,\lambda) = a_1^2+a_2^2+...a_n^2+2\rho_{12}a_1a_2+2\rho_{13}a_1a_3+.... - \lambda(a_1^2+a_2^2+...a_n^2-1) \end{equation}という総和を作り，各\(a\)で偏微分すると
\begin{equation} {\partial f\over \partial a_k}=0 \quad 1\le k \le n \end{equation} から以下の式が\(n\)本導出される． \begin{eqnarray} 2a_1+2\rho_{12}a_2+2\rho_{13}a_3+.... - 2\lambda a_1 &=& 0\nonumber\\ 2\rho_{21}a_1+2a_2+2\rho_{23}a_3+.... - 2\lambda a_2 &=& 0\nonumber\\ .....&& \nonumber\\ 2\rho_{n1}a_1+2\rho_{n2}a_2 + .....2a_n - 2\lambda a_n &=& 0 \label{eqs} \end{eqnarray}行列の形に書き換えると
\begin{equation} \left( \begin{array}{ccccc} 1 & \rho_{12} & \rho_{13} & \ldots &\\ \rho_{21} & 1 & \rho_{23} & \ldots &\\ ... & ... & ... & ... & ...\\ \rho_{i1} \ldots & \rho_{ik} & ... & ... &\\ ... & ... & ... & ... & ... \end{array} \right) \left( \begin{array}{c} a_{1} \\ a_2 \\ ...\\ a_{k} \\ ... \end{array} \right)= \lambda \left( \begin{array}{c} a_1 \\ a_2 \\ ...\\ a_k\\ ... \end{array} \right) \end{equation}とまとめることができる．この行列を見ると，相関係数からなる「相関行列」の固有ベクトルが求めるaの組み合わせ，未定乗数λがその固有ベクトルの寄与ということになる．さらに，式\eqref{eqs}の第一行目の両辺に\(a_1\)，2行目の両辺に\(a_2\)，．．．第\(n\)行目に\(a_n\)を掛けて２で割ると
\begin{eqnarray} a_1^2+\rho_{12}a_1a_2+\rho_{13}a_1a_3+.... &=& \lambda a_1^2 \nonumber\\ \rho_{21}a_2a_1+a_2^2+\rho_{23}a_2a_3+.... &=&\lambda a_2^2 \nonumber\\ .....&& \nonumber\\ \rho_{n1}a_na_1+\rho_{n2}a_na_2 + .....a_n^2 &=& \lambda a_n^2 \end{eqnarray}これらの式をすべて足し合わせると，式（8）と（9）から
\begin{equation} S = \lambda \end{equation}となることがわかる．よって，それぞれの主成分に対応して生じる分散は固有値となることがわかる．したがって固有値が大きければ大きいほどデータのばらつきの説明に寄与していることになる．一方で，全ての固有値の和は今の場合 \(n\) （固有値の総和と行列の対角和は等しい）であり，\(n\)が全ての固有ベクトルに対する分散の和となる．したがって，\(\lambda/n\) が，その固有値の寄与率となる．
大きい順に2つの固有値と対応する固有ベクトルを求め，それぞれの主成分の寄与率を求めよ．また，それぞれの主成分に対応する固有ベクトルの値を観察することで，この集団の学力がどのような科目の組み合わせで特徴づけられているか議論せよ．なお，固有ベクトルは長さ1に規格化されていることに注意．
宣言文 z(n,m) x(n,m), xsigma(n), xav(n) 成績データの読み込み read(1,*) ( x(n,m), n = 1, 4 ) 各科目nについて，平均値xav(n)と標準偏差xsigma(n)を求める 式(2) 平均値と標準偏差を用いて，x(n,m)からz(n,m)を作る 式(1), (4) ZTZ/mを求め，相関行列の各要素を出力する 式(5) サブルーチンjacobiを用い，相関行列の固有値と対応する固有ベクトルを求める 固有値の大きさ/4が寄与率，固有ベクトルの成分が対応する主成分の要素を表す