[Scipy-svn] r3274 - in trunk/scipy/sandbox/models: . tests

scipy-svn@scip... scipy-svn@scip...
Tue Aug 28 14:24:06 CDT 2007


Author: jonathan.taylor
Date: 2007-08-28 14:23:59 -0500 (Tue, 28 Aug 2007)
New Revision: 3274

Added:
   trunk/scipy/sandbox/models/tests/test_bspline.py
Modified:
   trunk/scipy/sandbox/models/bspline.py
   trunk/scipy/sandbox/models/bspline_module.py
   trunk/scipy/sandbox/models/gam.py
   trunk/scipy/sandbox/models/glm.py
   trunk/scipy/sandbox/models/info.py
   trunk/scipy/sandbox/models/setup.py
   trunk/scipy/sandbox/models/smoothers.py
Log:
added proper docstrings to BSpline and SmoothingSpline classes


Modified: trunk/scipy/sandbox/models/bspline.py
===================================================================
--- trunk/scipy/sandbox/models/bspline.py	2007-08-28 16:54:29 UTC (rev 3273)
+++ trunk/scipy/sandbox/models/bspline.py	2007-08-28 19:23:59 UTC (rev 3274)
@@ -5,9 +5,51 @@
 from scipy.optimize import golden
 from scipy.sandbox.models import _bspline
 
+def _band2array(a, lower=0, symmetric=False, hermitian=False):
+    """
+    Take an upper or lower triangular banded matrix and return a
+    numpy array.
+
+    INPUTS:
+       a         -- a matrix in upper or lower triangular banded matrix
+       lower     -- is the matrix upper or lower triangular?
+       symmetric -- if True, return the original result plus its transpose
+       hermitian -- if True (and symmetric False), return the original
+                    result plus its conjugate transposed
+
+    """
+
+    n = a.shape[1]
+    r = a.shape[0]
+    _a = 0
+
+    if not lower:
+        for j in range(r):
+            _b = N.diag(a[r-1-j],k=j)[j:(n+j),j:(n+j)]
+            _a += _b 
+            if symmetric and j > 0: _a += _b.T
+            elif hermitian and j > 0: _a += _b.conjugate().T
+    else:
+        for j in range(r):
+            _b = N.diag(a[j],k=j)[0:n,0:n]
+            _a += _b 
+            if symmetric and j > 0: _a += _b.T
+            elif hermitian and j > 0: _a += _b.conjugate().T
+        _a = _a.T
+
+    return _a
+
+
 def _upper2lower(ub):
     """
     Convert upper triangular banded matrix to lower banded form.
+
+    INPUTS:
+       ub  -- an upper triangular banded matrix
+
+    OUTPUTS: lb
+       lb  -- a lower triangular banded matrix with same entries
+              as ub
     """
 
     lb = N.zeros(ub.shape, ub.dtype)
@@ -19,7 +61,14 @@
 
 def _lower2upper(lb):
     """
-    Convert upper triangular banded matrix to lower banded form.
+    Convert lower triangular banded matrix to upper banded form.
+
+    INPUTS:
+       lb  -- a lower triangular banded matrix
+
+    OUTPUTS: ub
+       ub  -- an upper triangular banded matrix with same entries
+              as lb
     """
 
     ub = N.zeros(lb.shape, lb.dtype)
@@ -31,8 +80,22 @@
 
 def _triangle2unit(tb, lower=0):
     """
-    Take a banded triangular matrix and return its diagonal and the unit matrix:
-    the banded triangular matrix with 1's on the diagonal.
+    Take a banded triangular matrix and return its diagonal and the
+    unit matrix: the banded triangular matrix with 1's on the diagonal,
+    i.e. each row is divided by the corresponding entry on the diagonal.
+
+    INPUTS:
+       tb    -- a lower triangular banded matrix
+       lower -- if True, then tb is assumed to be lower triangular banded,
+                in which case return value is also lower triangular banded.
+
+    OUTPUTS: d, b
+       d     -- diagonal entries of tb
+       b     -- unit matrix: if lower is False, b is upper triangular
+                banded and its rows of have been divided by d,
+                else lower is True, b is lower triangular banded
+                and its columns have been divieed by d.
+
     """
     
     if lower: d = tb[0].copy()
@@ -43,9 +106,19 @@
         l = _upper2lower(tb)
         return d, _lower2upper(l / d)
     
-def _trace_symbanded(a,b, lower=0):
+def _trace_symbanded(a, b, lower=0):
     """
-    Compute the trace(a*b) for two upper or lower banded real symmetric matrices.
+    Compute the trace(ab) for two upper or banded real symmetric matrices
+    stored either in either upper or lower form.
+
+    INPUTS:
+       a, b    -- two banded real symmetric matrices (either lower or upper)
+       lower   -- if True, a and b are assumed to be the lower half
+
+
+    OUTPUTS: trace
+       trace   -- trace(ab)
+    
     """
 
     if lower:
@@ -58,7 +131,12 @@
 
 def _zero_triband(a, lower=0):
     """
-    Zero out unnecessary elements of a real symmetric banded matrix.
+    Explicitly zero out unused elements of a real symmetric banded matrix.
+
+    INPUTS:
+       a   -- a real symmetric banded matrix (either upper or lower hald)
+       lower   -- if True, a is assumed to be the lower half
+       
     """
 
     nrow, ncol = a.shape
@@ -68,18 +146,36 @@
         for i in range(nrow): a[i,0:i] = 0.
     return a
 
-def _zerofunc(x):
-    return N.zeros(x.shape, N.float)
 
+class BSpline(object):
 
-class BSpline:
+    '''
 
-    """
-    knots should be sorted, knots[0] is lower boundary, knots[1] is upper boundary
-    knots[1:-1] are internal knots
-    """
+    Bsplines of a given order and specified knots.
 
-    def __init__(self, knots, order=4, coef=None, M=None, eps=0.0):
+    Implementation is based on description in Chapter 5 of
+
+    Hastie, Tibshirani and Friedman (2001). "The Elements of Statistical
+    Learning." Springer-Verlag. 536 pages.
+
+
+    INPUTS:
+       knots  -- a sorted array of knots with knots[0] the lower boundary,
+                 knots[1] the upper boundary and knots[1:-1] the internal
+                 knots.
+       order  -- order of the Bspline, default is 4 which yields cubic
+                 splines
+       M      -- number of additional boundary knots, if None it defaults
+                 to order
+       coef   -- an optional array of real-valued coefficients for the Bspline
+                 of shape (knots.shape + 2 * (M - 1) - order,).
+       x      -- an optional set of x values at which to evaluate the
+                 Bspline to avoid extra evaluation in the __call__ method
+
+    '''
+
+    def __init__(self, knots, order=4, M=None, coef=None, x=None):
+
         knots = N.squeeze(N.unique(N.asarray(knots)))
 
         if knots.ndim != 1:
@@ -89,10 +185,9 @@
         if M is None:
             M = self.m
         self.M = M
-#         if self.M < self.m:
-#             raise 'multiplicity of knots, M, must be at least equal to order, m'
 
-        self.tau = N.hstack([[knots[0]-eps]*(self.M-1), knots, [knots[-1]+eps]*(self.M-1)])
+        self.tau = N.hstack([[knots[0]]*(self.M-1), knots, [knots[-1]]*(self.M-1)])
+
         self.K = knots.shape[0] - 2
         if coef is None:
             self.coef = N.zeros((self.K + 2 * self.M - self.m), N.float64)
@@ -100,12 +195,64 @@
             self.coef = N.squeeze(coef)
             if self.coef.shape != (self.K + 2 * self.M - self.m):
                 raise ValueError, 'coefficients of Bspline have incorrect shape'
+        if x is not None:
+            self.x = x
 
-    def __call__(self, x):
-        b = N.asarray(self.basis(x)).T
+    def _setx(self, x):
+        self._x = x
+        self._basisx = self.basis(self._x)
+
+    def _getx(self):
+        return self._x
+
+    x = property(_getx, _setx)
+
+    def __call__(self, *args):
+        """
+        Evaluate the BSpline at a given point, yielding
+        a matrix B and return
+
+        B * self.coef
+
+
+        INPUTS:
+           args -- optional arguments. If None, it returns self._basisx,
+                   the BSpline evaluated at the x values passed in __init__.
+                   Otherwise, return the BSpline evaluated at the
+                   first argument args[0].
+                   
+        OUTPUTS: y
+           y    -- value of Bspline at specified x values
+
+        BUGS:
+           If self has no attribute x, an exception will be raised
+           because self has no attribute _basisx.
+
+        """
+
+        if not args:
+            b = self._basisx.T
+        else:
+            x = args[0]
+            b = N.asarray(self.basis(x)).T
         return N.squeeze(N.dot(b, self.coef)) 
-    
+        
     def basis_element(self, x, i, d=0):
+        """
+        Evaluate a particular basis element of the BSpline,
+        or its derivative.
+
+        INPUTS:
+           x  -- x values at which to evaluate the basis element
+           i  -- which element of the BSpline to return
+           d  -- the order of derivative
+
+        OUTPUTS: y
+           y  -- value of d-th derivative of the i-th basis element
+                 of the BSpline at specified x values
+           
+        """
+        
         x = N.asarray(x, N.float64)
         _shape = x.shape
         if _shape == ():
@@ -122,7 +269,26 @@
         v.shape = _shape
         return v
 
-    def basis(self, x, d=0, upper=None, lower=None):
+    def basis(self, x, d=0, lower=None, upper=None):
+        """
+        Evaluate the basis of the BSpline or its derivative.
+        If lower or upper is specified, then only
+        the [lower:upper] elements of the basis are returned.
+
+        INPUTS:
+           x     -- x values at which to evaluate the basis element
+           i     -- which element of the BSpline to return
+           d     -- the order of derivative
+           lower -- optional lower limit of the set of basis
+                    elements 
+           upper -- optional upper limit of the set of basis
+                    elements 
+ 
+        OUTPUTS: y
+           y  -- value of d-th derivative of the basis elements
+                 of the BSpline at specified x values
+           
+        """
         x = N.asarray(x)
         _shape = x.shape
         if _shape == ():
@@ -154,9 +320,38 @@
             v[-1] = N.where(N.equal(x, self.tau[-1]), 1, v[-1])
         return v
 
-    def gram(self, d=0, full=False):
+    def gram(self, d=0):
         """
-        Compute Gram inner product matrix.
+        Compute Gram inner product matrix, storing it in lower
+        triangular banded form.
+        
+        The (i,j) entry is
+
+        G_ij = integral b_i^(d) b_j^(d)  
+
+        where b_i are the basis elements of the BSpline and (d) is the
+        d-th derivative.
+
+        If d is a matrix then, it is assumed to specify a differential
+        operator as follows: the first row represents the order of derivative
+        with the second row the coefficient corresponding to that order.
+
+        For instance:
+
+        [[2, 3],
+         [3, 1]]
+
+        represents 3 * f^(2) + 1 * f^(3).
+        
+        INPUTS:
+           d    -- which derivative to apply to each basis element,
+                   if d is a matrix, it is assumed to specify
+                   a differential operator as above
+
+        OUTPUTS: gram
+           gram -- the matrix of inner products of (derivatives)
+                   of the BSpline elements
+
         """
         
         d = N.squeeze(d)
@@ -184,18 +379,54 @@
     method = "target_df"
     target_df = 5
     default_pen = 1.0e-03
+    optimize = True
 
-    def smooth(self, y, x=None, weights=None):
-        if self.method == "target_df":
-            self.fit_target_df(y, x=x, weights=weights, df=self.target_df)
-        elif self.method == "optimize_gcv":
-            self.fit_optimize_gcv(y, x=x, weights=weights)
+    '''
+    A smoothing spline, which can be used to smooth scatterplots, i.e.
+    a list of (x,y) tuples.
 
+    See fit method for more information.
+    
+    '''
+    
     def fit(self, y, x=None, weights=None, pen=0.):
+        """
+        Fit the smoothing spline to a set of (x,y) pairs.
+
+        INPUTS:
+           y       -- response variable
+           x       -- if None, uses self.x
+           weights -- optional array of weights 
+           pen     -- constant in front of Gram matrix
+
+        OUTPUTS: None
+           The smoothing spline is determined by self.coef,
+           subsequent calls of __call__ will be the smoothing spline.
+           
+        ALGORITHM:
+           Formally, this solves a minimization:
+
+           fhat = ARGMIN_f SUM_i=1^n (y_i-f(x_i))^2 + pen * int f^(2)^2
+
+           See Chapter 5 of 
+  
+           Hastie, Tibshirani and Friedman (2001). "The Elements of Statistical
+           Learning." Springer-Verlag. 536 pages.
+
+           for more details.
+
+        TODO:
+           Should add arbitrary derivative penalty instead of just
+           second derivative.
+        """
+
         banded = True
 
         if x is None:
-            x = self.tau[(self.M-1):-(self.M-1)] # internal knots
+            x = self._x
+            bt = self._basisx.copy()
+        else:
+            bt = self.basis(x)
 
         if pen == 0.: # can't use cholesky for singular matrices
             banded = False
@@ -204,7 +435,6 @@
             raise ValueError, 'x and y shape do not agree, by default x are \
                the Bspline\'s internal knots'
         
-        bt = self.basis(x)
         if pen >= self.penmax:
             pen = self.penmax
 
@@ -225,13 +455,16 @@
         y = y[mask]
 
         self.df_total = y.shape[0]
-        bty = N.dot(bt, _w * y)
+
+        bty = N.squeeze(N.dot(bt, _w * y))
         self.N = y.shape[0]
+
         if not banded:
             self.btb = N.dot(bt, bt.T)
-            _g = band2array(self.g, lower=1, symmetric=True)
+            _g = _band2array(self.g, lower=1, symmetric=True)
             self.coef, _, self.rank = L.lstsq(self.btb + pen*_g, bty)[0:3]
             self.rank = min(self.rank, self.btb.shape[0])
+            del(_g) 
         else:
             self.btb = N.zeros(self.g.shape, N.float64)
             nband, nbasis = self.g.shape
@@ -239,7 +472,8 @@
                 for k in range(min(nband, nbasis-i)):
                     self.btb[k,i] = (bt[i] * bt[i+k]).sum()
 
-                bty.shape = (1,bty.shape[0])
+            bty.shape = (1,bty.shape[0])
+            self.pen = pen
             self.chol, self.coef = solveh_banded(self.btb + 
 						 pen*self.g,
 						 bty, lower=1)
@@ -248,9 +482,24 @@
         self.resid = y * self.weights - N.dot(self.coef, bt)	    
         self.pen = pen
 
+        del(bty); del(mask); del(bt)
+
+    def smooth(self, y, x=None, weights=None):
+        
+        if self.method == "target_df":
+            if hasattr(self, 'pen'):
+                self.fit(y, x=x, weights=weights, pen=self.pen)
+            else:
+                self.fit_target_df(y, x=x, weights=weights, df=self.target_df)
+        elif self.method == "optimize_gcv":
+            self.fit_optimize_gcv(y, x=x, weights=weights)
+
+
     def gcv(self):
         """
         Generalized cross-validation score of current fit.
+
+        TODO: addin a reference to Wahba, and whoever else I used.
         """
 
         norm_resid = (self.resid**2).sum()
@@ -258,6 +507,8 @@
 
     def df_resid(self):
         """
+        Residual degrees of freedom in the fit.
+        
         self.N - self.trace()
 
         where self.N is the number of observations of last fit.
@@ -267,15 +518,18 @@
 
     def df_fit(self):
         """
-        = self.trace()
+        How many degrees of freedom used in the fit? 
 
-        How many degrees of freedom used in the fit? 
+        self.trace()
+
         """
         return self.trace()
 
     def trace(self):
         """
         Trace of the smoothing matrix S(pen)
+
+        TODO: addin a reference to Wahba, and whoever else I used.
         """
 
         if self.pen > 0:
@@ -285,20 +539,36 @@
         else:
             return self.rank
 
-    def fit_target_df(self, y, x=None, df=None, weights=None, tol=1.0e-03):
+    def fit_target_df(self, y, x=None, df=None, weights=None, tol=1.0e-03,
+                      apen=0, bpen=1.0e-03):
+
         """
         Fit smoothing spline with approximately df degrees of freedom
         used in the fit, i.e. so that self.trace() is approximately df.
 
+        Uses binary search strategy.
+
         In general, df must be greater than the dimension of the null space
         of the Gram inner product. For cubic smoothing splines, this means
         that df > 2.
 
+        INPUTS:
+           y       -- response variable
+           x       -- if None, uses self.x
+           df      -- target degrees of freedom
+           weights -- optional array of weights 
+           tol     -- (relative) tolerance for convergence 
+           apen    -- lower bound of penalty for binary search
+           bpen    -- upper bound of penalty for binary search
+
+        OUTPUTS: None
+           The smoothing spline is determined by self.coef,
+           subsequent calls of __call__ will be the smoothing spline.
+
         """
 
         df = df or self.target_df
 
-        apen, bpen = 0, 1.0e-03
         olddf = y.shape[0] - self.m
 
         if hasattr(self, "pen"):
@@ -312,6 +582,7 @@
                 apen, bpen = 0., self.pen
 
         while True:
+
             curpen = 0.5 * (apen + bpen)
             self.fit(y, x=x, weights=weights, pen=curpen)
             curdf = self.trace()
@@ -326,7 +597,7 @@
                 break
 
     def fit_optimize_gcv(self, y, x=None, weights=None, tol=1.0e-03,
-		     bracket=(0,1.0e-03)):
+                         brack=(-100,20)):
         """
 	Fit smoothing spline trying to optimize GCV.
 
@@ -336,6 +607,18 @@
 	It is probably best to use target_df instead, as it is 
 	sometimes difficult to find a bracketing interval.
 
+        INPUTS:
+           y       -- response variable
+           x       -- if None, uses self.x
+           df      -- target degrees of freedom
+           weights -- optional array of weights 
+           tol     -- (relative) tolerance for convergence 
+           brack   -- an initial guess at the bracketing interval
+
+        OUTPUTS: None
+           The smoothing spline is determined by self.coef,
+           subsequent calls of __call__ will be the smoothing spline.
+
         """
 
         def _gcv(pen, y, x):
@@ -343,32 +626,7 @@
             a = self.gcv()
             return a
 
-        a = golden(_gcv, args=(y,x), brack=(-100,20), tol=tol)
+        a = golden(_gcv, args=(y,x), brack=bracket, tol=tol)
 
 
-def band2array(a, lower=0, symmetric=False, hermitian=False):
-    """
-    Take an upper or lower triangular banded matrix and return a matrix using
-    LAPACK storage convention. For testing banded Cholesky decomposition, etc.
-    """
 
-    n = a.shape[1]
-    r = a.shape[0]
-    _a = 0
-
-    if not lower:
-        for j in range(r):
-            _b = N.diag(a[r-1-j],k=j)[j:(n+j),j:(n+j)]
-            _a += _b 
-            if symmetric and j > 0: _a += _b.T
-            elif hermitian and j > 0: _a += _b.conjugate().T
-    else:
-        for j in range(r):
-            _b = N.diag(a[j],k=j)[0:n,0:n]
-            _a += _b 
-            if symmetric and j > 0: _a += _b.T
-            elif hermitian and j > 0: _a += _b.conjugate().T
-        _a = _a.T
-
-    return _a
-

Modified: trunk/scipy/sandbox/models/bspline_module.py
===================================================================
--- trunk/scipy/sandbox/models/bspline_module.py	2007-08-28 16:54:29 UTC (rev 3273)
+++ trunk/scipy/sandbox/models/bspline_module.py	2007-08-28 19:23:59 UTC (rev 3274)
@@ -132,11 +132,11 @@
     PyArrayObject *basis;
     double *data;
 
-
     basis = (PyArrayObject *) PyArray_SimpleNew(2, dim, PyArray_DOUBLE);
     data = (double *) basis->data;
     bspline(&data, x, Nx[0], knots, Nknots[0], m, d, lower, upper);
     return_val = (PyObject *) basis;
+    Py_DECREF((PyObject *) basis); 
 
     '''    
 
@@ -184,7 +184,9 @@
     
     double bspline_quad(double *knots, int nknots,
                         int m, int l, int r, int dl, int dr) 
+
         /* This is based on scipy.integrate.fixed_quad */
+
     {
         double *y;
         double qx[%(nq)d]={%(qx)s};
@@ -202,8 +204,8 @@
         lower = l - m - 1; 
 	if (lower < 0) { lower = 0;}
 	upper = lower + 2 * m + 4;
-	if (upper > nknots - 1) {upper = nknots-1;}
-/*	upper = nknots - m; */
+	if (upper > nknots - 1) { upper = nknots-1; }
+
         for (k=lower; k<upper; k++) {
             partial = 0.;
             a = knots[k]; b=knots[k+1];
@@ -216,7 +218,10 @@
             for (kk=0; kk<nq; kk++) {
                 partial += y[kk] * qw[kk];
             }
+            free(y); /* bspline_prod malloc's memory, but does not free it */
+
             result += (b - a) * partial / 2.;
+
         }
 
         return(result);
@@ -260,6 +265,7 @@
     data = (double *) gram->data;
     bspline_gram(&data, knots, Nknots[0], m, dl, dr);
     return_val = (PyObject *) gram;
+    Py_DECREF((PyObject *) gram); 
 
     '''    
 
@@ -324,7 +330,8 @@
     data = (double *) invband->data;
     invband_compute(&data, L, NL[1], NL[0]-1);
 
-     return_val = (PyObject *) invband;
+    return_val = (PyObject *) invband;
+    Py_DECREF((PyObject *) invband); 
 
     '''    
 

Modified: trunk/scipy/sandbox/models/gam.py
===================================================================
--- trunk/scipy/sandbox/models/gam.py	2007-08-28 16:54:29 UTC (rev 3273)
+++ trunk/scipy/sandbox/models/gam.py	2007-08-28 19:23:59 UTC (rev 3274)
@@ -1,6 +1,8 @@
 """
 Generalized additive models
+
 """
+
 import numpy as N
 
 from scipy.sandbox.models import family
@@ -12,7 +14,7 @@
     _x.sort()
     n = x.shape[0]
     # taken form smooth.spline in R
-    print "herenow"
+
     if n < 50:
         nknots = n
     else:
@@ -29,7 +31,8 @@
         else:
             nknots = 200 + (n - 3200.)**0.2
         knots = _x[N.linspace(0, n-1, nknots).astype(N.int32)]
-    s = SmoothingSpline(knots)
+
+    s = SmoothingSpline(knots, x=x.copy())
     s.gram(d=2)
     s.target_df = 5
     return s
@@ -41,7 +44,7 @@
         self.offset = offset
 
     def __call__(self, *args, **kw):
-        return self.fn(*args, **kw) + offset
+        return self.fn(*args, **kw) + self.offset
 
 class results:
 
@@ -62,7 +65,7 @@
         return N.sum(self.smoothed(design), axis=0) + self.alpha
 
     def smoothed(self, design):
-        return N.array([self.smoothers[i](design[:,i]) + self.offset[i] for i in range(design.shape[1])])
+        return N.array([self.smoothers[i]() + self.offset[i] for i in range(design.shape[1])])
 
 class additive_model:
 
@@ -89,16 +92,16 @@
         offset = N.zeros(self.design.shape[1], N.float64)
         alpha = (Y * self.weights).sum() / self.weights.sum()
         for i in range(self.design.shape[1]):
-            tmp = self.smoothers[i](self.design[:,i])
-            self.smoothers[i].smooth(Y - alpha - mu + tmp, x=self.design[:,i],
+            tmp = self.smoothers[i]()
+            self.smoothers[i].smooth(Y - alpha - mu + tmp,
 				     weights=self.weights)
-            tmp2 = self.smoothers[i](self.design[:,i])
+            tmp2 = self.smoothers[i]()
             offset[i] = -(tmp2*self.weights).sum() / self.weights.sum()
             mu += tmp2 - tmp
 
         return results(Y, alpha, self.design, self.smoothers, self.family, offset)
 
-    def cont(self, tol=1.0e-02):
+    def cont(self, tol=1.0e-04):
 	
         curdev = (((self.results.Y - self.results.predict(self.design))**2) * self.weights).sum()
 
@@ -124,9 +127,9 @@
         offset = N.zeros(self.design.shape[1], N.float64)
 
         for i in range(self.design.shape[1]):
-            self.smoothers[i].smooth(Y - alpha - mu, x=self.design[:,i],
+            self.smoothers[i].smooth(Y - alpha - mu,
 				     weights=self.weights)
-            tmp = self.smoothers[i](self.design[:,i])
+            tmp = self.smoothers[i]()
             offset[i] = (tmp * self.weights).sum() / self.weights.sum()
             tmp -= tmp.sum()
             mu += tmp
@@ -190,8 +193,7 @@
         return self.results
 
 
-if __name__ == "__main__":
-
+def _run():
     import numpy.random as R
     n = lambda x: (x - x.mean()) / x.std()
     n_ = lambda x: (x - x.mean()) 
@@ -220,11 +222,11 @@
     toc = time.time()
     m.fit(b)
     tic = time.time()
-    import pylab
-    pylab.figure(num=1)
-    pylab.plot(x1, n(m.smoothers[0](x1))); pylab.plot(x1, n(f1(x1)), linewidth=2)
-    pylab.figure(num=2)
-    pylab.plot(x2, n(m.smoothers[1](x2))); pylab.plot(x2, n(f2(x2)), linewidth=2); 
+##     import pylab
+##     pylab.figure(num=1)
+##     pylab.plot(x1, n(m.smoothers[0](x1)), 'r'); pylab.plot(x1, n(f1(x1)), linewidth=2)
+##     pylab.figure(num=2)
+##     pylab.plot(x2, n(m.smoothers[1](x2)), 'r'); pylab.plot(x2, n(f2(x2)), linewidth=2); 
     print tic-toc
 
     f = family.Poisson()
@@ -235,9 +237,12 @@
     m.fit(p)
     tic = time.time()
     print tic-toc
-    pylab.figure(num=1)
-    pylab.plot(x1, n(m.smoothers[0](x1))); pylab.plot(x1, n(f1(x1)), linewidth=2) 
-    pylab.figure(num=2)
-    pylab.plot(x2, n(m.smoothers[1](x2))); pylab.plot(x2, n(f2(x2)), linewidth=2)
-    pylab.show()
+##     pylab.figure(num=1)
+##     pylab.plot(x1, n(m.smoothers[0](x1)), 'b'); pylab.plot(x1, n(f1(x1)), linewidth=2) 
+##     pylab.figure(num=2)
+##     pylab.plot(x2, n(m.smoothers[1](x2)), 'b'); pylab.plot(x2, n(f2(x2)), linewidth=2)
+##     pylab.show()
 
+
+if __name__ == "__main__":
+    _run()

Modified: trunk/scipy/sandbox/models/glm.py
===================================================================
--- trunk/scipy/sandbox/models/glm.py	2007-08-28 16:54:29 UTC (rev 3273)
+++ trunk/scipy/sandbox/models/glm.py	2007-08-28 19:23:59 UTC (rev 3274)
@@ -1,5 +1,6 @@
 """
-General linear model
+General linear models
+--------------------
 """
 import numpy as N
 from scipy.sandbox.models import family
@@ -14,7 +15,7 @@
         self.weights = 1
         self.initialize(design)
 
-    def __iter__(self):
+    def __iter__(self):                                                   
         self.iter = 0
         self.dev = N.inf
         return self

Modified: trunk/scipy/sandbox/models/info.py
===================================================================
--- trunk/scipy/sandbox/models/info.py	2007-08-28 16:54:29 UTC (rev 3273)
+++ trunk/scipy/sandbox/models/info.py	2007-08-28 19:23:59 UTC (rev 3274)
@@ -6,7 +6,7 @@
 
   - `ols_model` (ordinary least square regression)
   - `wls_model` (weighted least square regression)
-  - `ar_model` (autoregression)
+  - `ar_model` (autoregressive model)
 
  - `glm.model` (generalized linear models)
  - robust statistical models

Modified: trunk/scipy/sandbox/models/setup.py
===================================================================
--- trunk/scipy/sandbox/models/setup.py	2007-08-28 16:54:29 UTC (rev 3273)
+++ trunk/scipy/sandbox/models/setup.py	2007-08-28 19:23:59 UTC (rev 3274)
@@ -8,6 +8,8 @@
     config.add_data_dir('tests')
 
     try:
+        import sys
+        print sys.path
         from scipy.sandbox.models.bspline_module import mod
         n, s, d = weave_ext(mod)
         config.add_extension(n, s, **d)

Modified: trunk/scipy/sandbox/models/smoothers.py
===================================================================
--- trunk/scipy/sandbox/models/smoothers.py	2007-08-28 16:54:29 UTC (rev 3273)
+++ trunk/scipy/sandbox/models/smoothers.py	2007-08-28 19:23:59 UTC (rev 3274)
@@ -10,7 +10,7 @@
 from scipy.optimize import golden
 
 from scipy.sandbox.models import _bspline
-from scipy.sandbox.models.bspline import bspline, band2array
+from scipy.sandbox.models.bspline import bspline, _band2array
 
 
 class poly_smoother:
@@ -103,7 +103,7 @@
         self.N = y.shape[0]
         if not banded:
             self.btb = N.dot(bt, bt.T)
-            _g = band2array(self.g, lower=1, symmetric=True)
+            _g = _band2array(self.g, lower=1, symmetric=True)
             self.coef, _, self.rank = L.lstsq(self.btb + pen*_g, bty)[0:3]
             self.rank = min(self.rank, self.btb.shape[0])
         else:

Added: trunk/scipy/sandbox/models/tests/test_bspline.py
===================================================================
--- trunk/scipy/sandbox/models/tests/test_bspline.py	2007-08-28 16:54:29 UTC (rev 3273)
+++ trunk/scipy/sandbox/models/tests/test_bspline.py	2007-08-28 19:23:59 UTC (rev 3274)
@@ -0,0 +1,23 @@
+"""
+Test functions for models.glm
+"""
+
+import numpy as N
+from numpy.testing import NumpyTest, NumpyTestCase
+
+import scipy.sandbox.models as S
+import scipy.sandbox.models.bspline as B
+
+
+class test_BSpline(NumpyTestCase):
+
+    def test1(self):
+        b = B.BSpline(N.linspace(0,10,11), x=N.linspace(0,10,101))
+        old = b._basisx.shape
+        b.x = N.linspace(0,10,51)
+        new = b._basisx.shape
+        self.assertEqual((old[0], 51), new)
+
+
+if __name__ == "__main__":
+    NumpyTest().run()



More information about the Scipy-svn mailing list