Juliaで主成分分析

周りの人と読んでいる本(ベイズ推論による機械学習)の影響でJuliaをさわってみました。
まだあまりむずかしいことができないので、操作に慣れるために主成分分析をしてみました。
（主成分分析を一発で実行する関数もあるかもしれませんが、練習なので教科書通りやってみます）

using PyCall, PyPlot
using Distributions

# データの作成
n = 300
μ　 = [1.0, 3.0]
Σ　　= [10 3
        3 2]
X = rand(MvNormal( μ, Σ), n);

MvNormal関数で2次元のガウス分布に従うデータを作成しました。 全角の文字も変数にできるのでギリシャ文字を変数にできて楽しい。

# データから平均を計算して、データから引く
mu = [mean(X[1, :])
            mean(X[2, :])]

X_m = copy(X)
X_m[1, :] = X[1, :] - mu[1];
X_m[2, :] = X[2, :] - mu[2];

2x300のデータから2x1の形のベクトルをうまく引き算する方法がわからなかったのでダサい方法で。。

# 共分散行列の計算
Σ = 1.0 / n *  X_m * X_m'

2×2 Array{Float64,2}: 10.5301 3.24329 3.24329 2.1232

# 固有値分解
λ, S = eig(Σ)
println(λ)
println(S)

# 線形変換
y = S' * X_m;

[1.01741, 11.6359] [0.322704 -0.9465; -0.9465 -0.322704]

# 可視化用の第一成分、第二成分
p_ax1 = [mu[1] mu[1] + S[1,1]*4
                mu[2] mu[2] + S[2,1] *4]
p_ax2 = [mu[1] mu[1] + S[1,2]*4
                mu[2] mu[2] + S[2,2]*4 ]

2×2 Array{Float64,2}: 1.15528 -2.63072 3.10447 1.81366

figure("result")
clf()
subplot(121)
title("PCA result")
plot(X[1,:], X[2, :], "+")
plot(p_ax1[1, :], p_ax1[2, :])
plot(p_ax2[1, :], p_ax2[2, :])

axis("equal")

subplot(122)
title("Linear Transformed")
plot(y[1,:], y[2, :], "+")
axis("equal")

show()

プロットもPyPlotを使えば（pythonになじみがあれば）あまり迷うことなくできました。 次は画像を使って見ようかな、データ用意するの面倒だけど。