[SciPy-User] orthonormal polynomials (cont.)

[Charles R Harris]

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.
>
>
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.

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