因子分析の＜沼＞に片足だけ突っ込んでみた５

『マンガでわかる統計学　因子分析編』で解説される手順をエクセルExcelで確かめながら行う、因子分析の学習の続きです。前回は、Step12まで進めました。 

Step13 元の行列の共通性の値を、Step12で求めた因子負荷量の積の行列における対角要素の値に置き換える

Step14 Step13における行列から固有値と固有ベクトルを求める

Step15 任意の共通性の値が１を超えそうになるまで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 "２５回反復しましたが、収束しませんでした。", vbExclamation
        NumberOfTimes = 0
    Else
        Call yacobi4FA
    End If
    
End Sub

改変ポイントを簡単に解説すると、

Subプロシージャに変えました。それに伴い、引数を定数にしています。手抜きですが、継続的に使うものでもないので。orz
定数nの値は、行列の次元数に合わせて適宜変更してください。

因子負荷量とその積を求める数式をセルに入力しています。

因子負荷量の積の計算後の共通性の値が１を超えたかどうか判定しています。

共通性の値が１を超えていなければ、固有ベクトルの値を別のセルに書き出します。

元の行列に新しい共通性の値をコピペして、再帰でループを回します。
２５回反復しても収束しないときはループを止めるようにしてありますが、反復回数は適宜変更してください。

で、Step10で求めた行列をコピペし直して、上のマクロを実行した結果がこちらです。

*1:オリジナルのyacobi関数については以下をご参照ください。
エクセルで操る！ヤコビ法による固有値計算、固有ベクトル計算