因子分析の<沼>に片足だけ突っ込んでみた5
『マンガでわかる統計学 因子分析編』で解説される手順をエクセルExcelで確かめながら行う、因子分析の学習の続きです。前回は、Step12まで進めました。
Step13 元の行列の共通性の値を、Step12で求めた因子負荷量の積の行列における対角要素の値に置き換える
Step14 Step13における行列から固有値と固有ベクトルを求める
Step15 任意の共通性の値が1を超えそうになるまでStep12からStep14までを繰り返す
繰り返しということなので、作業を自動化すべくyacobi関数を少し弄ってみました。*1
Option Explicit
Dim NumberOfTimes As Integer
'ヤコビ法による固有値解析 For 因子分析
Sub yacobi4FA() 'n As Long)
Const n = 6
Dim i As Long
Dim j As Long
Dim k As Long
Dim im As Long
Dim jm As Long
Dim amax As Double
Dim Q As Double
Dim a() As Double
Dim sr() As Double
Dim Anchor As Boolean
ReDim a(n, n)
ReDim sr(n, n)
'初期設定
For i = 1 To n
For j = 1 To n
a(i, j) = Cells(i, j).Value
Next
Next
For i = 1 To n
For j = 1 To n
If i = j Then
sr(i, j) = 1
Else
sr(i, j) = 0
End If
Next
Next
For i = 1 To n
For j = 1 To n
Cells(i + 3 * n, j).Value = sr(i, j)
Next
Next
If NumberOfTimes = 0 Then
Cells(3, 2 + n).FormulaR1C1 = "=ATAN2(R[-2]C[-1]-R[-1]C,2*R[-2]C)/2"
Range(Cells(2 * n + 1, 1), Cells(3 * n, n)).FormulaArray = _
"=MMULT(MMULT(RC[" & n & "]:R[" & n - 1 & "]C[" & 2 * n - 1 & "],R[-" & n & _
"]C:R[-1]C[" & n - 1 & "]),R[-" & n & "]C[" & n & "]:R[-1]C[" & 2 * n - 1 & "])"
Range(Cells(2 * n + 1, n + 1), Cells(3 * n, 2 * n)).FormulaArray = _
"=TRANSPOSE(R[-" & n & "]C:R[-1]C[" & n - 1 & "])"
Range(Cells(3 * n + 1, n + 1), Cells(4 * n, 2 * n)).FormulaArray = _
"=MMULT(RC[-" & n & "]:R[" & n - 1 & "]C[-1],R[-" & 2 * n & "]C:R[-" & n + 1 & "]C[" & n - 1 & "])"
End If
' 最大値の探索
For k = 1 To 100
amax = 0#
For i = 1 To n - 1
For j = i + 1 To n
If Abs(a(i, j)) > amax Then
amax = Abs(a(i, j))
im = i
jm = j
End If
Next
Next
If amax < 0.0000001 Then Exit For
' Arの設定
For i = 1 To n
For j = 1 To n
Cells(i + n, j).Value = a(i, j)
Next
Next
' 最大値の取出し
Cells(1, 1 + n).Value = a(im, im)
Cells(1, 2 + n).Value = a(im, jm)
Cells(2, 1 + n).Value = a(jm, im)
Cells(2, 2 + n).Value = a(jm, jm)
' Srの設定
Q = Cells(3, 2 + n).Value
For i = 1 To n
For j = 1 To n
If i = j Then
sr(i, j) = 1
Else
sr(i, j) = 0
End If
Next
Next
sr(im, im) = Cos(Q)
sr(im, jm) = -Sin(Q)
sr(jm, im) = Sin(Q)
sr(jm, jm) = Cos(Q)
For i = 1 To n
For j = 1 To n
Cells(i + n, j + n).Value = sr(i, j)
Next
Next
' Sのコピー
For i = 1 To n
For j = 1 To n
sr(i, j) = Cells(i + 3 * n, j + n).Value
Next
Next
For i = 1 To n
For j = 1 To n
Cells(i + 3 * n, j).Value = sr(i, j)
Next
Next
' Ar+1の取出し
For i = 1 To n
For j = 1 To n
a(i, j) = Cells(i + 2 * n, j).Value
Next
Next
Next
If NumberOfTimes = 0 Then
For i = 1 To n
Cells(4 * n + i, 1).FormulaR1C1 = "=SQRT(R[" & n * -2 - i + 1 & "]C1)*R[" & -n & "]C[" & n & "]"
Cells(4 * n + i, 2).FormulaR1C1 = "=SQRT(R[" & n * -2 - i + 2 & "]C2)*R[" & -n & "]C[" & n & "]"
Next
Range(Cells(4 * n + 1, 3), Cells(5 * n, n + 2)).FormulaArray = _
"=MMULT(RC[-2]:R[" & n - 1 & "]C[-1],TRANSPOSE(RC1:R[" & n - 1 & "]C2))"
End If
For i = 1 To n
If Cells(4 * n + 1, 3).Offset(i - 1, i - 1).Value > 1 Then
Anchor = True
Exit For
End If
Next
If Anchor = True Then Exit Sub
Range(Cells(5 * n + 2, 1), Cells(6 * n + 1, 2)).Value = _
Range(Cells(4 * n + 1, 1), Cells(5 * n + 1, 2)).Value
For i = 1 To n
Cells(i, i).Value = Cells(4 * n + 1, 3).Offset(i - 1, i - 1).Value
Next
NumberOfTimes = NumberOfTimes + 1
If NumberOfTimes = 25 Then
MsgBox "25回反復しましたが、収束しませんでした。", vbExclamation
NumberOfTimes = 0
Else
Call yacobi4FA
End If
End Sub
改変ポイントを簡単に解説すると、
Subプロシージャに変えました。それに伴い、引数を定数にしています。手抜きですが、継続的に使うものでもないので。orz
定数nの値は、行列の次元数に合わせて適宜変更してください。
因子負荷量とその積を求める数式をセルに入力しています。
因子負荷量の積の計算後の共通性の値が1を超えたかどうか判定しています。
共通性の値が1を超えていなければ、固有ベクトルの値を別のセルに書き出します。
元の行列に新しい共通性の値をコピペして、再帰でループを回します。
25回反復しても収束しないときはループを止めるようにしてありますが、反復回数は適宜変更してください。
で、Step10で求めた行列をコピペし直して、上のマクロを実行した結果がこちらです。
*1:オリジナルのyacobi関数については以下をご参照ください。
エクセルで操る!ヤコビ法による固有値計算、固有ベクトル計算