[Scipy-svn] r3274 - in trunk/scipy/sandbox/models: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Tue Aug 28 15:24:06 EDT 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