[SciPy-User] orthonormal polynomials (cont.)

[Charles R Harris]

On Tue, Jun 14, 2011 at 3:57 PM, <josef.pktd at gmail.com> wrote:

> On Tue, Jun 14, 2011 at 5:48 PM, Charles R Harris
> <charlesr.harris at gmail.com> wrote:
> >
> >
> > On Tue, Jun 14, 2011 at 3:15 PM, Charles R Harris
> > <charlesr.harris at gmail.com> wrote:
> >>
> >>
> >> On Tue, Jun 14, 2011 at 3:01 PM, <josef.pktd at gmail.com> wrote:
> >>>
> >>> On Tue, Jun 14, 2011 at 4:40 PM, Charles R Harris
> >>> <charlesr.harris at gmail.com> wrote:
> >>> >
> >>> >
> >>> > On Tue, Jun 14, 2011 at 12:10 PM, nicky van foreest
> >>> > <vanforeest at gmail.com>
> >>> > wrote:
> >>> >>
> >>> >> Hi,
> >>> >>
> >>> >> Without understanding the details... I recall from numerical recipes
> >>> >> in C that Gram Schmidt is a very risky recipe. I don't know whether
> >>> >> this advice also pertains to fitting polynomials, however,
> >>>
> >>> I read some warnings about the stability of Gram Schmidt, but the idea
> >>> is that if I can choose the right weight function, then we should need
> >>> only a few polynomials. So, I guess in that case numerical stability
> >>> wouldn't be very relevant.
> >>>
> >>> >>
> >>> >> Nicky
> >>> >>
> >>> >> On 14 June 2011 18:58,  <josef.pktd at gmail.com> wrote:
> >>> >> > (I'm continuing the story with orthogonal polynomial density
> >>> >> > estimation, and found a nice new paper
> >>> >> >
> http://www.informaworld.com/smpp/content~db=all~content=a933669464 )
> >>> >> >
> >>> >> > Last time I managed to get orthonormal polynomials out of scipy
> with
> >>> >> > weight 1, and it worked well for density estimation.
> >>> >> >
> >>> >> > Now, I would like to construct my own orthonormal polynomials for
> >>> >> > arbitrary weights. (The weights represent a base density around
> >>> >> > which
> >>> >> > we make the polynomial expansion).
> >>> >> >
> >>> >> > The reference refers to Gram-Schmidt or Emmerson recurrence.
> >>> >> >
> >>> >> > Is there a reasonably easy way to get the polynomial coefficients
> >>> >> > for
> >>> >> > this with numscipython?
> >>> >> >
> >>> >
> >>> > What do you mean by 'polynomial'? If you want the values of a set of
> >>> > polynomials orthonormal on a given set of points, you want the 'q' in
> a
> >>> > qr
> >>> > factorization of a (row) weighted Vandermonde matrix.  However, I
> would
> >>> > suggest using a weighted chebvander instead for numerical stability.
> >>>
> >>> Following your suggestion last time to use QR, I had figured out how
> >>> to get the orthonormal basis for a given set of points.
> >>> Now, I would like to get the functional version (not just for a given
> >>> set of points), that is an orthonormal polynomial basis like Hermite,
> >>> Legendre, Laguerre and Jacobi, only for any kind of weight function,
> >>> where the weight function is chosen depending on the data.
> >>>
> >>
> >> But in what basis? The columns of the inverse of 'R' in QR will give you
> >> the orthonormal polynomials as series in whatever basis you used for the
> >> columns of the pseudo-Vandermonde matrix.
>
> with basis you mean hear for example the power series  (x**i,
> i=0,1,..)  that's what they use, but there is also a reference to
> using fourier polynomials which I haven't looked at for this case.
>
> >>
> >
> > Example.
> >
> > In [1]: from numpy.polynomial.polynomial import polyvander
> >
> > In [2]: v = polyvander(linspace(-1, 1, 1000), 3)
> >
> > In [3]: around(inv(qr(v, mode='r'))*sqrt(1000./2), 5)
> > Out[3]:
> > array([[-0.70711,  0.     ,  0.79057, -0.     ],
> >        [ 0.     ,  1.22352, -0.     , -2.80345],
> >        [ 0.     ,  0.     , -2.36697,  0.     ],
> >        [ 0.     ,  0.     ,  0.     ,  4.66309]])
> >
> > The columns are approx. the coefficients of the normalized Legendre
> > functions as a power series up to a sign.
>
> Looks interesting. It will take me a while to figure out what this
> does, but I think I get the idea.
>
> Thanks,
>
>
The normalization factor comes from integrating the columns of q, i.e., \int
p^2 ~= dt*\sum_i (q_{ij})^2 = 2/1000. I really should have weighted the
first and last rows of the Vandermonde matrix by 1/sqrt(2) so that the
integration was trapazoidal, but you get the idea.

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20110614/8bb28b31/attachment.html>