射影追跡をニュートン法で実装した

本当にメモ程度です。

射影追跡とは

　標本の非ガウス方向を見つけるアルゴリズムのこと（これが、混合した音声などの独立成分分析に応用できる）

射影追跡の考え方

１、ガウス分布の尖度は3である。→尖度が3から遠いほど、非ガウス
２、勾配法で尖度を更新してもよいが、ニュートン法の方が速い
３、外れ値に強くするため、非ガウス性尺度を定義する

射影追跡の手順

１、標本を中心化&球状化（白色化）
２、非ガウス方向のパラメータが収束まで更新

実装

GitHubにコードをあげています。
また、中心化と球状化をしなくてもいいように標本を作成しています。

import numpy as np

# 標本数
N = 1000

# Gauss分布を生成
x_gauss = np.random.randn(N)

# 一様分布を生成
x_uniform = np.random.rand(N) * 5.0 - 2.5

x = np.array([x_gauss, x_uniform])


def nomalize(b):
    return b / np.sqrt(b[0]**2 + b[1]**2)


def main():
    # g(s) = s^3 のとき
    # 初期のbを仮定
    b = nomalize(np.array([1, 0]))
    threshold = 0.0000000001
    diff = 1

    while diff > threshold:
        print(b)
        tmp_1, tmp_2 = 0, 0
        for i in range(N):
            tmp_1 += 3 * (np.dot(b, x[:, i]) ** 2)
            tmp_2 += x[:, i] * (np.dot(b, x[:, i]) ** 3)

        new_b = tmp_2 / N - tmp_1 * b / N
        diff = (b[0] - nomalize(new_b)[0])**2 + (b[1] - nomalize(new_b)[1])**2
        b = nomalize(new_b)


if __name__ == '__main__':
    main()

結果

確かに非ガウス方向に射影できている。