[SciPy-User] B-spline basis functions?
Charles R Harris
Wed Aug 1 12:21:59 CDT 2012
On Wed, Aug 1, 2012 at 10:53 AM, Nathaniel Smith <firstname.lastname@example.org> wrote:
> On Tue, Jul 31, 2012 at 4:46 PM, Charles R Harris
> <email@example.com> wrote:
> > On Tue, Jul 31, 2012 at 9:37 AM, Nathaniel Smith <firstname.lastname@example.org> wrote:
> >> Hi all,
> >> I'd like to be able to do spline regression in patsy, which means
> >> that I need to be able to compute b-spline basis functions. I am not
> >> an initiate into the mysteries of practical spline computations, but I
> >> *think* the stuff in scipy.signal is not quite usable as is, because
> >> it's focused on doing interpolation directly rather than exposing the
> >> basis functions themselves?
> >> Specifically, to achieve feature parity with R , I need to be able to
> >> take
> >> - an arbitrary order
> >> - an arbitrary collection of knot positions (which may be irregularly
> >> spaced)
> >> - a vector x of points at which to evaluate the basis functions
> >> and spit out the value of each spline basis function evaluated at each
> >> point in the x vector.
> >> It looks like scipy.signal.bspline *might* be useful, but I can't
> >> quite tell? Or alternatively someone might have some code lying around
> >> to do this already?
> >> Basically I have a copy of Schumaker here and I'm hoping someone will
> >> save me from having to read it :-).
> > I have this floating around
> > v[i] = spl.splev(x, (knots, d[i], deg))
> This looks fabulous, thank you! From a quick look it seems to be
> producing numerically identical results to R's spline.des().
> The bit with the coefficient vectors like [0, 0, 0, 1, 0] worries me a
> bit -- do you know if it's producing the entire spline basis
> internally on every call, and then throwing away most of them when
> forming the inner product? If so then in practice it's probably not a
> show-stopper, but it would be a pointless factor-of-len(knots)
> increase in running time for no reason.
Yes it is. It was a quick and dirty routine and producing the array wasn't
the bottleneck. Also, since the b-splines themselves have small support, a
lot of extra zero values are produced. It wouldn't be difficult to improve
things quite a bit using a hacked version of the splev.f routine (De Boor's
algorithm <http://en.wikipedia.org/wiki/De_Boor%27s_algorithm>). IIRC, the
b-splines in signal are periodic b-splines with reflected boundary
conditions, so aren't as general. I could be wrong about that, however.
-------------- next part --------------
An HTML attachment was scrubbed...
More information about the SciPy-User