FastICA

fastica

NB. Fast ICA
NB. http://en.wikipedia.org/wiki/FastICA
require 'math/lapack math/lapack/geev'
require 'plot'
center =: -"1 +/%# 			NB. subtract column means
mp     =: +/ .*			NB. matrix product
cov    =: (mp"2 1)~&.|: ] % <:^:(1&<)@#	NB. Covariance (note, does NOT center)
I      =: =@:i.			NB. unit matrix
clean  =: * |@* 			NB. clean near-to-zero values
whiten =: (mp ([ mp (I@# % %:)@] mp |:@[)&>/@(2b110&geev_jlapack_)@cov)@center
NB. single component extraction
f1  =: ^.@(6&o.)
g1  =: 7&o.	
g1p =: 1 - *:@(7&o.)
f2  =: -@^@-@-:@*:
g2  =: * ^@-@-:@*:
g2p =: (1-*:) * ^@-@-:@*:

f=:f1
g=:g1
gp=:g1p
NB. inputs:
NB. C: number of components to extract
NB. X: N x M data N samples, M dimensions, C < M!
NB. outputs:
NB. W: Unmixing matrix, here each row projects X onto into independent component.
NB. R: Independent components matrix, with M columns representing a sample with C dimensions.
fastica =: 4 : 0
'M N'=. $y
w =. (x,N)$0
for_p. i.x do.
  wp   =. (% +/) ? N $ 0 
  wold =.   N $ 0
  ct =. 0
  while. (ct < 10000) *. -. wold -: wp do.
    wold =. wp
    wp =. - M %~ ((g y mp wp) mp y) - wp * +/ (gp y mp wp)
    if. p>0 do.
      wp =. wp - +/ (wp mp w{~p-1)
    end.
    wp =. (% +/&.:*:) wp
    ct =. ct+1
  end.
  w =. wp p} w
end.
w;y mp |:w
)

NB. Example: separate sine from noise:

n   =: center 0.1 * ? 1000 3 $ 0
s1  =:   1 o. 10%~i.&.(10&*) 100
s2  =: * 1 o.3%~i.&.(10&*) 100
s3  =: center 50 %~ 100 | i. 1000

mix =: (%"1 +/) ? 3 3 $0 
in  =: n + ( s1 ,. s2 ,. s3) mp mix
win  =: whiten in
f   =: 3  fastica win
pd'reset'
pd'multi  1 3'
pd'ygroup 0 0 0'
pd ('key recov1 recov2 recov3';|: 1{:: f) ; ('key mix1 mix2 mix3'; |: in) ,&< ('key orig1 orig2 orig3'; s1 , s2 ,: s3)
pd'show'