covariance_iir.py 2017-01-14 Jan

covariance_iir.py

#!/usr/bin/env python2

from __future__ import division
import numpy as np

__version__ = "2017-01-13 jan  denis-bz-py t-online de"

#...............................................................................
class Covariance_iir( object ):
    """ running Covariance_iir filter, up-weighting more recent data like IIR

    Example:
        covariance_filter = Covariance_iir( a=0.1  or halflife=7 )
        for x in data:  # e.g. 3d
            cov = covariance_filter( x )  # 3 x 3

    https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance

    FIR the first `nfir` data points, then IIR with factor `a`: see iir.py
    """

    def __init__( s, halflife=7, a=None, nfir=1 ):
        if a is None:
            a = 0.7 / halflife  # a, 1 - a ?
        assert 0 < a < 1, a
        s.a = a
        if nfir is None:
            nfir = 0.7 / a
        assert nfir > 0, nfir
        s.nfir = nfir
        s.t = -1

    def __call__( s, x ):
        """ state `av += a * (x - av)`, fraction `a` of the way towards `x`
            `deltadelta = (x - old av) * (x - new av)`, * = outer product
            state `cov += a * (deltadelta - cov)`
        """
        s.t += 1
        if s.t == 0:
            s.av = np.copy( x )  # init x0, not 0
            s.cov = np.zeros( (len(x), len(x)), dtype=x.dtype )
        else:
                # a_t = 1 1/2 1/3 ... FIR coefs, then IIR a --
            a_t = 1. / (1 + s.t) if s.t < s.nfir \
                else s.a
            delta = x - s.av
            s.av += a_t * delta
            delta2 = x - s.av
            s.cov *= (1 - a_t)
            s.cov += a_t * np.outer( delta, delta2 )

        return np.copy( s.cov )

    def __str__( s ):
        return "Covariance_iir-halflife%.3g-nfir%.0f" % (.7 / s.a, s.nfir)

normalcloud.py

#!/usr/bin/env python2
from __future__ import division
import numpy as np

def normalcloud( N, scale=[1,1], tilt=30, add=[0,0] ):
    """ -> N x 2 cloud scale tilt add """
    scale = np.asarray(scale)
    add = np.asarray(add)
    cloud = np.random.multivariate_normal( [0,0], np.diag(scale**2), (N,) )
    return np.dot( cloud, rotate(np.radians(tilt)) )  + add

def rotate(rad):
    """ rotate . [1 0] -> [cos sin] """
    c, s = np.cos(rad), np.sin(rad)
    return np.array( [[c, -s],
                    [s, c]] )

def normalcloud2( nx=1, ny=1 ):
    """ 2 normal clouds in unit square, ~ L """
    nxy = np.array([nx,ny])
    n3 = int( N * 1.3 / 3 )
    clouds = np.vstack((
        normalcloud( n3, nxy * np.r_[ .25, .05 ], 10, [.5, .1] ),
        normalcloud( 2*n3, nxy * np.r_[ .4, .1 ], 80, [.3, .6] )))
    clouds01 = clouds[ np.all( (0 <= clouds) & (clouds <= 1), axis=1 )]
    # print len(clouds01)
    return clouds01[:N]

#...............................................................................
if __name__ == "__main__":
    import sys
    import pylab as pl
    from etc import numpyutil as nu, pca  # ~/py/bz/etc/pca.py

    np.set_printoptions( threshold=20, edgeitems=10, linewidth=140,
            formatter = dict( float = lambda x: "%.4g" % x ))  # float arrays %.4g

    N = 100
    tilt = 0
    plot = 0
    seed = 0

    # to change these params, run this.py a=1 b=None 'c = ...'  in sh or ipython
    for arg in sys.argv[1:]:
        exec( arg )
    np.random.seed(seed)
    print "\nN %d  tilt %g" % (N, tilt)

    cloud = normalcloud( N, scale=[10,100], tilt=tilt )  # N, 2
    print "cov:\n", np.cov(cloud.T)
    x, y = cloud.T
    jx = x.argsort()
    print "x:", nu.ints(x[jx])
    print "y:", nu.ints(y[jx])
    xlen = nu.av2(x)
    ylen = nu.av2(y)
    print "|y| %.0f / |x| %.0f = %.2g" % (ylen, xlen, ylen / xlen)
        # tilt 0:  |y| 102 / |x| 10 
        # tilt 30: |y| 89 / |x| 52 = 1.7
    print ""

    p = pca.PCA( cloud )
    print "pca Vt:\n", p.Vt
    projx, projy = p.project( cloud ).T
    print "projx:", nu.minavmax( projx )
    print "projy:", nu.minavmax( projy )

    if plot:
        pl.figure( figsize=[10,10] )
        pl.plot( projx, projy, "+" )
        pl.show()

test-covar.py

# test-covar.py
#!/usr/bin/env python2

from __future__ import division
import sys
import numpy as np

from etc.filter.covariance_iir import Covariance_iir
from normalcloud import normalcloud
    # def normalcloud( N, scale=[1,1], tilt=30, add=[0,0] ):

np.set_printoptions( threshold=20, edgeitems=10, linewidth=140,
        formatter = dict( float = lambda x: "%.4g" % x ))  # float arrays %.4g

N = 100
tilt = 30
halflife = 7
nfir = np.inf
plot = 0
seed = 0

    # to change these params, run this.py a=1 b=None 'c = ...'  in sh or ipython
for arg in sys.argv[1:]:
    exec( arg )
np.random.seed(seed)
print "\ncloud: N %d  tilt %g" % (N, tilt)

#...............................................................................
covfilter = Covariance_iir( nfir=nfir, halflife=halflife )

cloud = normalcloud( N, scale=[10,100], tilt=tilt )  # N, 2
for x in cloud:
    cov = covfilter( x )
print "%s: \n%s" % (covfilter, cov )
print "np.cov:\n", np.cov( cloud.T, bias=True )  # bias: / N not N-1