Go言語の数値計算ライブラリGonumを使う(2)特異値分解（SVD）、主成分分析（PCA)を計算してみる。

今回の事例はこちら。

Visual C# (C_sharp)の数値計算ライブラリ MathNET Numericsを使う(8)　特異値分解（SVD）、主成分分析（PCA)を計算してみる（ちょうど奥村先生が記事を出されてたので 

元データはこちら。

https://okumuralab.org/~okumura/stat/biplot.html

GonumはSvdも普通に計算できる。主成分分析そのものはないが、Svdから計算できる。

コードはこちら。

package main

import (

    "encoding/csv"

    "fmt"

    "log"

    "math"

    "os"

    "strconv"

    "gonum.org/v1/gonum/mat"

)

func main() {

    file, err := os.Open("atest2014chu.csv")

    if err != nil {

        log.Fatal(err)

    }

    defer file.Close()

    r := csv.NewReader(file)

    rows, err := r.ReadAll()

    if err != nil {

        log.Fatal(err)

    }

    prefectures := 47

    ntest := 4

    data := make([]float64, prefectures*ntest)

    count := 0

    for i := 1; i < prefectures+1; i++ {

        for j := 1; j < ntest+1; j++ {

            data[count], err = strconv.ParseFloat(rows[i][j], 64)

            count++

        }

    }

    m := mat.NewDense(prefectures, ntest, data)

    // fmt.Println(mat.Formatted(m))

    colmean := make([]float64, ntest)

    for j := 0; j < ntest; j++ {

        colmean[j] = 0.0

        for i := 0; i < prefectures; i++ {

            colmean[j] += m.At(i, j)

        }

        colmean[j] /= float64(prefectures)

    }

    for i := 0; i < prefectures; i++ {

        for j := 0; j < ntest; j++ {

            m.Set(i, j, m.At(i, j)-colmean[j])

        }

    }

    var svd mat.SVD

    ok := svd.Factorize(m, mat.SVDFull)

    if !ok {

        panic("SVD factorization failed")

    }

    var s []float64

    s = svd.Values(s)

    fmt.Printf("Singular values: %6.4v\n", s)

    scale := 1.0 / math.Sqrt(float64(prefectures-1))

    fmt.Printf("Principle Components:  %6.4v, %6.4v, %6.4v, %6.4v\n", s[0]*scale, s[1]*scale, s[2]*scale, s[3]*scale)

    fmt.Println()

    u := mat.NewDense(prefectures, prefectures, nil)

    vt := mat.NewDense(ntest, ntest, nil)

    uw := mat.NewDense(prefectures, ntest, nil)

    w := mat.NewDense(prefectures, ntest, nil)

    w.Set(0, 0, s[0])

    w.Set(1, 1, s[1])

    w.Set(2, 2, s[2])

    w.Set(3, 3, s[3])

    svd.UTo(u)

    svd.VTo(vt)

    uw.Mul(u, w)

    fmt.Println("Vt = ")

    fmt.Printf("%6.4v\n", mat.Formatted(vt))

    fmt.Println("U * w=")

    fmt.Printf("%6.4v\n", mat.Formatted(uw))

}

実行結果。奥村先生の結果とも合っている。

Singular values: [ 30.96 7.878 4.667 2.474]
Principle Components: 4.565, 1.162, 0.6881, 0.3647

Vt =
⎡ 0.3236 0.4347 -0.5039 -0.6726⎤
⎢ 0.4297 0.6693 0.0467 0.6044⎥
⎢ 0.5571 -0.5789 -0.5225 0.2854⎥
⎣ 0.6326 -0.1671 0.6862 -0.3177⎦
U * w=
⎡ -1.492 0.02216 0.5907 -0.8963⎤
⎢ 2.589 -0.004574 -0.9497 -0.1746⎥
⎢ -2.58 3.189 -0.1859 -0.4263⎥
⎢ -0.5212 2.05 0.4443 -0.3466⎥
⎢ 10.42 1.072 -1.125 -0.6339⎥
⎢ 1.239 0.9562 -0.1983 -0.1714⎥
⎢ -3.603 1.611 -0.4861 -0.07854⎥
⎢ -1.033 1.824 -0.9083 0.869⎥
⎢ -0.4417 0.4782 0.1723 0.3661⎥
⎢ 4.713 1.252 0.4587 0.4556⎥
⎢ -0.7564 0.9939 0.4923 0.1595⎥
⎢ 0.2437 0.8785 0.5878 -0.1001⎥
⎢ 3.425 0.7737 0.2736 0.2602⎥
⎢ 0.5736 0.1931 1.204 0.04571⎥
⎢ 0.02382 -0.8335 -0.1593 -0.3068⎥
⎢ 7.179 1.109 -0.4105 0.4169⎥
⎢ 6.21 0.1469 -0.3071 -0.1538⎥
⎢ 11.62 -0.4558 -0.1334 0.2942⎥
⎢ 0.1286 1.291 0.2788 0.04524⎥
⎢ -1.764 -0.6595 -1.103 -0.6448⎥
⎢ 2.859 -0.916 -0.1945 0.5845⎥
⎢ 5.529 -1.185 0.3971 -0.2346⎥
⎢ 3.153 -1.654 0.9602 0.276⎥
⎢ -2.415 -1.642 -0.07539 0.1647⎥
⎢ -4.133 -1.138 -1.093 0.5753⎥
⎢ 0.4995 -0.1421 0.2482 0.2125⎥
⎢ -5.568 -1.831 0.4811 -0.4051⎥
⎢ 2.393 -1.359 -0.1824 -0.08378⎥
⎢ 0.854 -1.181 0.2829 0.2805⎥
⎢ -4.871 -1.961 -0.2454 -0.2042⎥
⎢ 0.06871 0.07501 0.07541 0.1263⎥
⎢ -0.4648 1.768 0.02756 -0.0413⎥
⎢ -5.202 -0.7719 -0.9768 -0.2363⎥
⎢ 1.197 -0.5776 -0.2145 -0.4274⎥
⎢ 4.423 -0.7638 0.001467 -0.04943⎥
⎢ -0.997 -1.318 -1.562 0.2369⎥
⎢ 0.5026 -0.1598 0.1824 0.4005⎥
⎢ 3.764 0.2062 1.264 -0.09054⎥
⎢ -8.657 0.4732 -0.6591 0.1549⎥
⎢ -3.18 -0.1144 0.1917 -0.01086⎥
⎢ -5.367 0.3199 0.2999 -0.2689⎥
⎢ -1.29 -0.7063 0.1801 -0.3674⎥
⎢ 1.152 -0.393 1.636 -0.08071⎥
⎢ -2.165 0.3837 -1.283 -0.1773⎥
⎢ -0.03332 -1.685 0.1978 0.4835⎥
⎢ -3.168 -0.6218 0.7835 -0.3305⎥
⎣ -15.06 1.007 0.7406 0.5332⎦

なかなか使えるな。次はFFTかな。