ScipyでICA

#!/usr/bin/env python
# -*- coding:utf-8 -*-
import scipy as sp
import numpy as np
import matplotlib.pyplot as plt

def g( val, alpha = 0.01 ):
    return np.tanh( alpha * val )

def gd( val, alpha = 0.01 ):
    return alpha * ( 1.0 - np.asarray( np.tanh( alpha * val ) ) )

def sym_orth( sig ):
    a,b,c = np.linalg.svd( sig * sig.T )
    return np.dot( np.dot( a, np.dot( np.diag( np.sqrt( 1./ b ) ), c ) ), sig )

def fast_ica( sig ):
    m,n = sig.shape

    #白色化
    u,s,v = np.linalg.svd( sig, full_matrices=False)
    v = np.sqrt( n ) * v

    #分離フィルタの初期化
    w = np.random.rand( m, m )
    w = w / np.matrix( [ [ np.linalg.norm( w[ :, x ] ) ] * m for x in xrange( m ) ] ).T
    w = sym_orth( w )

    #ICA
    err = 1.0
    it = 0
    max_iter = 3000
    while 1e-12 < min( err , it < max_iter ):
        last_w = w
        w1 =  np.dot( v , g( np.dot( w.T , v ) ).T )
        w2 = np.multiply( np.dot( gd( np.dot( w.T, v ) ), np.ones( (n, 1 ) ) ) , w )
        w = ( w1 - w2 ) / n
        w = sym_orth( w )
        err = np.sum( np.abs( np.abs(w) - np.abs(last_w) ) )
        it += 1

    #0番目の信号のゲインを基準に分離フィルタを作成
    W = np.dot( np.diag( np.ravel( u[0,] ) ), np.dot( np.diag( s ) , w.T ) )
    return np.dot( W, v )  / np.sqrt( n )

def main():
    #スーパーガウシアンな分布だしラプラス分布で
    s = np.random.laplace( 0., 1.0, ( 2, 10000 ) )
    a = np.asmatrix( np.random.rand( 2,2 ) )
    x = a * s
    y = np.asarray( fast_ica( x ) )
    plt.scatter( y[0,], y[1,] )
    plt.show()
    return

if __name__ == '__main__':
    main()