From 902806c1d8b8184f7c811f119edaa67117b24ca3 Mon Sep 17 00:00:00 2001 From: flatberg Date: Wed, 31 Oct 2007 18:19:00 +0000 Subject: [PATCH] Sandbox version of models --- pyblm/crossvalidation.py | 13 +- pyblm/models.py | 282 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 291 insertions(+), 4 deletions(-) create mode 100644 pyblm/models.py diff --git a/pyblm/crossvalidation.py b/pyblm/crossvalidation.py index a824068..ab7aded 100644 --- a/pyblm/crossvalidation.py +++ b/pyblm/crossvalidation.py @@ -140,7 +140,7 @@ def pca_jk(a, aopt, n_blocks=None): def pls_jk(X, Y, a_opt, nsets=None, center=True, verbose=False): """ Returns jack-knife segements of W. - *Parameters*: + *Parameters*: X : {array} Main data matrix (m, n) @@ -251,12 +251,14 @@ def find_aopt_from_sep(sep, method='vanilla'): closely before deciding on the optimal number of components. *Parameters*: + sep : {array} Squared error of prediction method : ['vanilla', '75perc'] Mehtod used to estimate optimal number of components *Returns*: + aopt : {integer} A guess on the optimal number of components """ @@ -297,12 +299,14 @@ def cv(N, K, randomise=True, sequential=False): Use sequential sampling *Returns*: + training : {array-like} training-indices validation : {array-like} validation-indices *Notes*: + If randomise is true, a copy of index is shuffled before partitioning, otherwise its order is preserved in training and validation. @@ -334,9 +338,10 @@ def cv(N, K, randomise=True, sequential=False): def diag_cv(shape, nsets=9): """Generates K (training, validation) index pairs. - Parameters: + *Parameters*: + N : {integer} - alpha -- scalar, approx. portion of data perturbed + alpha -- scalar, approx. portion of data perturbed """ try: m, n = shape @@ -371,7 +376,7 @@ def class_error(Yhat, Y, method='vanilla'): return Yhat_c, err -def class_errorII(T, Y, method='lda'): +def _class_errorII(T, Y, method='lda'): """ Not used ... """ pass diff --git a/pyblm/models.py b/pyblm/models.py new file mode 100644 index 0000000..17c6730 --- /dev/null +++ b/pyblm/models.py @@ -0,0 +1,282 @@ +"""Bilinear models""" + +from numpy import expand_dims + +from engines import pca + +def mean_center(x, axis=0): + """Mean center across axis.""" + return expand_dims(-x.mean(axis), axis) + +def scale(x, axis=0): + """ Scale across axis.""" + scale = ones((x.shape[axis-1],)) + #scale = 1./x.std(axis) + return expand_dims(scale, axis) + +class Model(object): + def __init__(name="johndoe"): + self.name = name + self.options = {} + + def save(self, filename='pca.ml'): + pass + + def load(self): + pass + + def clear(self): + for param in self.__dict__.keys(): + if param.startswith("_") and param[1]!="_": + exec "del self." + param + + def clear_core(self): + for param in self.__dict__.keys(): + if param.startswith("_"): + exec "del self." + param + #self.clear() + + +class PCA(Model): + def __init__(self, x, amax=10): + Model.__init__(self, name="PCA") + self._x = x + self.amax = amax + self.aopt = amax + + # properties + def amax(): + doc = "maximum number of components" + def fget(self): + return self._amax + def fset(self, a): + assert(a>0) + a = min(a, min(self.x.shape)) + if hasattr(self, "_amax"): + if a>self._amax: + # fixme: do a model update + raise NotImplementedError + self._amax = a + def fdel(self): + pass + return locals() + amax = property(**amax()) + + def tot_var(): + doc = "total variance" + def fget(self): + if not hasattr(self, "_tot_var"): + self._tot_var = (self.xw**2).sum() + return self._tot_var + def fset(self, tv): + self._tot_var = tv + def fdel(self): + del self._tot_var + return locals() + tot_var = property(**tot_var()) + + def scores(): + doc = "pca scores" + def fget(self): + if not hasattr(self, "_scores"): + u, s, v, tot_var = pcaengine(self.xw, self.amax) + self._scores = u + self.singvals = s + self.loadings = v + self.tot_var = tot_var + return self._scores[:,:self.amax] + def fset(self, t): + self._scores = t + def fdel(self): + del self._scores + return locals() # credit: David Niergarth + scores = property(**scores()) + + def loadings(): + doc = "pca loadings" + def fget(self): + if not hasattr(self, "_loadings"): + u, s, v, tot_var = pcaengine(self.xw, self.amax) + self._loadings = v + self.scores = u + self.singvals = s + self.tot_var = tot_var + return self._loadings[:,:self.amax] + def fdel(self): + del self._loadings + def fset(self, p): + self._loadings = p + return locals() + loadings = property(**loadings()) + + def singvals(): + doc = "Singular values" + def fget(self): + if not hasattr(self, "_singvals"): + u, s, v, tot_var = pcaengine(self.xw, self.amax) + self._singvals = s + self.scores = u + self.loadings = v + self.tot_var = tot_var + return self._singvals[:self.amax] + def fset(self, w): + self._singvals = w + def fdel(self): + del self._singvals + return locals() + singvals = property(**singvals()) + + def x(): + doc = "x is readonly, may not be deleted" + def fget(self): + return self._x + def fdel(self): + pass + return locals() + x = property(**x()) + + def xadd(): + doc = "column means" + def fget(self): + if not hasattr(self, "_xadd"): + self._xadd = center(self.x, axis=0) + return self._xadd + def fset(self, mnx): + if hasattr(self, "_xc"): + del self._xc + self._xadd = mnx + def fdel(self): + del self._xadd + if hasattr(self, "_xc"): + del self._xc + return locals() + xadd = property(**xadd()) + + def xc(): + doc = "centered input data" + def fget(self): + if not hasattr(self, "_xc"): + self._xc = self.x + self.xadd + return self._xc + def fset(self, xc): + self._xc = xc + def fdel(self): + del self._xc + return locals() + xc = property(**xc()) + + def xw(): + doc = "scaled input data" + def fget(self): + if not hasattr(self, "_xw"): + if self.x.shape==self.row_metric.shape: + self._xw = dot(self.xc, self.row_metric) + else: + self._xw = self.xc * self.row_metric + return self._xw + def fset(self, xw): + self._xw = xw + def fdel(self): + del self._xw + return locals() + xw = property(**xw()) + + def explained_variance(): + doc = "explained variance" + def fget(self): + if not hasattr(self, "_explained_variance"): + self._explained_variance = 100*(self.singvals**2)/self.tot_var + return self._explained_variance[:self.amax] + def fset(self, ev): + self._explained_variance = ev + def fdel(self): + del self._explained_variance + return locals() + explained_variance = property(**explained_variance()) + + def residuals(): + doc = "residuals" + def fget(self): + if not hasattr(self, "_residuals"): + res = empty((self.amax, self.x.shape[0], self.x.shape[1])) + for a in range(self.amax): + res[a,:,:] = self.xw - dot(self.scores[:,:a+1], self.loadings[:,:a+1].T) + self._residuals = res + return self._residuals + def fset(self, e): + self._residuals = e + def fdel(self): + del self._residuals + return locals() + residuals = property(**residuals()) + + def leverage(): + doc = "objects leverage" + def fget(self): + if not hasattr(self, "_leverage"): + u = self.scores/self.singvals + self._leverage = empty((self.amax, u.shape[0])) + for i in range(self.amax): + self._leverage[i,:] = 1./u.shape[0] + (u[:,:i+1]**2).sum(1) + return self._leverage[:self.amax,:] + def fset(self, lev): + self._leverage = lev + def fdel(self): + del self._leverage + return locals() + leverage = property(**leverage()) + + def row_metric(): + doc = "row metric" + def fget(self): + if not hasattr(self, "_row_metric"): + self._row_metric = scale(self.xc, axis=0) + return self._row_metric + def fset(self, w): + self._row_metric = w + def fdel(self): + del self._row_metric + if hasattr(self, "_xd"): + del self.xd + return locals() + row_metric = property(**row_metric()) + + def column_metric(): + doc = "column metric" + def fget(self): + if not hasattr(self, "_column_metric"): + self._column_metric = scale(self.xc, axis=1) + return self._column_metric + def fset(self, w): + + self._column_metric = w + # update model + def fdel(self): + del self._column_metric + if hasattr(self, "_xd"): + del self.xd + return locals() + column_metric = property(**column_metric()) + + def blm_update(self, a, b): + pass + + def append_columns(self, cols): + pass + + def append_rows(self, rows): + pass + + def delete_columns(self, index): + pass + + def delete_rows(self, index): + pass + + def reweight(self, ) + + +if __name__ == "__main__": + X = random.rand(4,10) + pcaobj = PCA(X) + print "explained variance" + str(pcaobj.explained_variance)