[SciPy-User] B-spline basis functions?

Charles R Harris charlesr.harris@gmail....
Wed Aug 1 12:21:59 CDT 2012


On Wed, Aug 1, 2012 at 10:53 AM, Nathaniel Smith <njs@pobox.com> wrote:

> On Tue, Jul 31, 2012 at 4:46 PM, Charles R Harris
> <charlesr.harris@gmail.com> wrote:
> > On Tue, Jul 31, 2012 at 9:37 AM, Nathaniel Smith <njs@pobox.com> wrote:
> >>
> >> Hi all,
> >>
> >> I'd like to be able to do spline regression in patsy[1], 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 [2], 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.

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-user/attachments/20120801/3375f35b/attachment.html 


More information about the SciPy-User mailing list