gistfile1.txt

NB.  Signal Processing in J
NB. ========================
NB. Mainly written as exercise, so no guarantees
NB. cocurrent <'jimproc'
NB. Assumptions : 1-cell = 1 sample of dimensionsionality $ samp
NB. set of samples in 2-cell

NB. Monad mean y
NB. Mean of data samples in y
mean =: +/ % #

NB. Helpers to cov
NB. n-1 if n> 1, else n
NB. nMinusOne =. ((]`<:)@.(>&1)) @ #
nMinusOne =: <:^:(1&<)@#

NB. Difference with sample mean
diffMean  =: (-"1 (+/ % #))

NB. Matrix Square
MatrixSquare =: +/ . *"2 1 ~

NB. Monad cov y
NB. Covariance of datasamples in y
cov  =:  (MatrixSquare &. |:@:diffMean % nMinusOne) 
NB.    sum over matrix multiplication of per sample differences with the mean over all samples, divided by n-1 if n>1, else n;
NB. x nancov y
NB.   x value given to "missing value"
NB.   y data
NB.   calculate an estimate for the covariance matrix
nancov =: cov@(] #~ -.@((+./"1)@:=))

NB. x pca y
NB.   x = value given to missing values
NB.   y = data to be transformed
NB.   returns a boxed array containing:
NB.   0) the coefficients (new basis in columns)
NB.   1) data projected onto this new basis
NB.   2) the variance explained, cumulative
NB.   3) the column means (needed when transforming back)
pca =: 4 : 0
require 'math/lapack'
require 'math/lapack/geev'
c=. x nancov y
'vec val bar' =. geev_jlapack_ c
NB. column means
cm =. x (+/%#)@(] #~ -.@((+./"1)@:=)) y
vec;(vec +/ .*~ y -"1 cm);((+/\%+/)val);cm
)

a =: ? 1000 10 $ 0
res =: _ pca a
coeff =: 0 {:: res
PC1to6 =: 1 {:: res
varExplained =: 2{::res 
colMeans =: 3{:: res NB. for transforming back