[SciPy-User] orthogonal polynomials ?

Mon May 16 22:16:20 EDT 2011

On Mon, May 16, 2011 at 7:47 PM, <josef.pktd at gmail.com> wrote:

> On Mon, May 16, 2011 at 9:29 PM, Charles R Harris
> <charlesr.harris at gmail.com> wrote:
> >
> >
> > On Mon, May 16, 2011 at 6:12 PM, <josef.pktd at gmail.com> wrote:
> >>
> >> On Mon, May 16, 2011 at 5:58 PM, Charles R Harris
> >> <charlesr.harris at gmail.com> wrote:
> >> >
> >> >
> >> > On Mon, May 16, 2011 at 12:40 PM, <josef.pktd at gmail.com> wrote:
> >> >>
> >> >> On Mon, May 16, 2011 at 12:47 PM,  <josef.pktd at gmail.com> wrote:
> >> >> > On Mon, May 16, 2011 at 12:29 PM, Charles R Harris
> >> >> > <charlesr.harris at gmail.com> wrote:
> >> >> >>
> >> >> >>
> >> >> >> On Mon, May 16, 2011 at 10:22 AM, <josef.pktd at gmail.com> wrote:
> >> >> >>>
> >> >> >>> On Mon, May 16, 2011 at 11:27 AM, Charles R Harris
> >> >> >>> <charlesr.harris at gmail.com> wrote:
> >> >> >>> >
> >> >> >>> >
> >> >> >>> > On Mon, May 16, 2011 at 9:17 AM, Charles R Harris
> >> >> >>> > <charlesr.harris at gmail.com> wrote:
> >> >> >>> >>
> >> >> >>> >>
> >> >> >>> >> On Mon, May 16, 2011 at 9:11 AM, <josef.pktd at gmail.com>
> wrote:
> >> >> >>> >>>
> >> >> >>> >>> On Sat, May 14, 2011 at 6:04 PM, Charles R Harris
> >> >> >>> >>> <charlesr.harris at gmail.com> wrote:
> >> >> >>> >>> >
> >> >> >>> >>> >
> >> >> >>> >>> > On Sat, May 14, 2011 at 2:26 PM, <josef.pktd at gmail.com>
> >> >> >>> >>> > wrote:
> >> >> >>> >>> >>
> >> >> >>> >>> >> On Sat, May 14, 2011 at 4:11 PM, nicky van foreest
> >> >> >>> >>> >> <vanforeest at gmail.com>
> >> >> >>> >>> >> wrote:
> >> >> >>> >>> >> > Hi,
> >> >> >>> >>> >> >
> >> >> >>> >>> >> > Might this be what you want:
> >> >> >>> >>> >> >
> >> >> >>> >>> >> > The first eleven probabilists' Hermite polynomials are:
> >> >> >>> >>> >> >
> >> >> >>> >>> >> > ...
> >> >> >>> >>> >> >
> >> >> >>> >>> >> > My chromium browser does not seem to paste pngs. Anyway,
> >> >> >>> >>> >> > check
> >> >> >>> >>> >> >
> >> >> >>> >>> >> >
> >> >> >>> >>> >> > http://en.wikipedia.org/wiki/Hermite_polynomials
> >> >> >>> >>> >> >
> >> >> >>> >>> >> > and you'll see that the first polynomial is 1, the
> second
> >> >> >>> >>> >> > x,
> >> >> >>> >>> >> > and
> >> >> >>> >>> >> > so
> >> >> >>> >>> >> > forth. From my courses on quantum mechanics I recall
> that
> >> >> >>> >>> >> > these
> >> >> >>> >>> >> > polynomials are, with respect to some weight function,
> >> >> >>> >>> >> > orthogonal.
> >> >> >>> >>> >>
> >> >> >>> >>> >> Thanks, I haven't looked at that yet, we should add
> >> >> >>> >>> >> wikipedia
> >> >> >>> >>> >> to
> >> >> >>> >>> >> the
> >> >> >>> >>> >> scipy.special docs.
> >> >> >>> >>> >>
> >> >> >>> >>> >> However, I would like to change the last part "with
> respect
> >> >> >>> >>> >> to
> >> >> >>> >>> >> some
> >> >> >>> >>> >> weight function"
> >> >> >>> >>> >>
> >> >> >>> >>> >>
> http://en.wikipedia.org/wiki/Hermite_polynomials#Orthogonality
> >> >> >>> >>> >>
> >> >> >>> >>> >> Instead of Gaussian weights I would like uniform weights
> on
> >> >> >>> >>> >> bounded
> >> >> >>> >>> >> support. And I have never seen anything about changing the
> >> >> >>> >>> >> weight
> >> >> >>> >>> >> function for the orthogonal basis of these kind of
> >> >> >>> >>> >> polynomials.
> >> >> >>> >>> >>
> >> >> >>> >>> >
> >> >> >>> >>> > In numpy 1.6, you can use the Legendre polynomials. They
> are
> >> >> >>> >>> > orthogonal
> >> >> >>> >>> > on
> >> >> >>> >>> > [-1,1] as has been mentioned, but can be mapped to other
> >> >> >>> >>> > domains.
> >> >> >>> >>> > For
> >> >> >>> >>> > example
> >> >> >>> >>> >
> >> >> >>> >>> > In [1]: from numpy.polynomial import Legendre as L
> >> >> >>> >>> >
> >> >> >>> >>> > In [2]: for i in range(5): plot(*L([0]*i + [1],
> >> >> >>> >>> > domain=[0,1]).linspace())
> >> >> >>> >>> >    ...:
> >> >> >>> >>> >
> >> >> >>> >>> > produces the attached plots.
> >> >> >>> >>>
> >> >> >>> >>> I'm still on numpy 1.5 so this will have to wait a bit.
> >> >> >>> >>>
> >> >> >>> >>> > <snip>
> >> >> >>> >>> >
> >> >> >>> >>> > Chuck
> >> >> >>> >>> >
> >> >> >>> >>>
> >> >> >>> >>>
> >> >> >>> >>> as a first application for orthogonal polynomials I was
> trying
> >> >> >>> >>> to
> >> >> >>> >>> get
> >> >> >>> >>> an estimate for a density, but I haven't figured out the
> >> >> >>> >>> weighting
> >> >> >>> >>> yet.
> >> >> >>> >>>
> >> >> >>> >>> Fourier polynomials work better for this.
> >> >> >>> >>>
> >> >> >>> >>
> >> >> >>> >> You might want to try Chebyshev then, the Cheybyshev
> >> >> >>> >> polynomialas
> >> >> >>> >> are
> >> >> >>> >> essentially cosines and will handle the ends better. Weighting
> >> >> >>> >> might
> >> >> >>> >> also
> >> >> >>> >> help, as I expect the distribution of the errors are somewhat
> >> >> >>> >> Poisson.
> >> >> >>> >>
> >> >> >>> >
> >> >> >>> > I should mention that all the polynomial fits will give you the
> >> >> >>> > same
> >> >> >>> > results, but the Chebyshev fits are more numerically stable.
> The
> >> >> >>> > general
> >> >> >>> > approach is to overfit, i.e., use more polynomials than needed
> >> >> >>> > and
> >> >> >>> > then
> >> >> >>> > truncate the series resulting in a faux min/max approximation.
> >> >> >>> > Unlike
> >> >> >>> > power
> >> >> >>> > series, the coefficients of the Cheybshev series will tend to
> >> >> >>> > decrease
> >> >> >>> > rapidly at some point.
> >> >> >>>
> >> >> >>> I think I might have still something wrong with the way I use the
> >> >> >>> scipy.special polynomials
> >> >> >>>
> >> >> >>> for a large sample size with 10000 observations, the nice graph
> is
> >> >> >>> fourier with 20 elements, the second (not so nice) is with
> >> >> >>> scipy.special.chebyt with 500 polynomials. The graph for 20
> >> >> >>> Chebychev
> >> >> >>> polynomials looks very similar
> >> >> >>>
> >> >> >>> Chebychev doesn't want to get low enough to adjust to the low
> part.
> >> >> >>>
> >> >> >>> (Note: I rescale to [0,1] for fourier, and [-1,1] for chebyt)
> >> >> >>>
> >> >> >>
> >> >> >> That certainly doesn't look right. Could you mail me the data
> >> >> >> offline?
> >> >> >> Also,
> >> >> >> what are you fitting, the histogram, the cdf, or...?
> >> >> >
> >> >> > The numbers are generated (by a function in scikits.statsmodels).
> >> >> > It's
> >> >> > all in the script that I posted, except I keep changing things.
> >> >> >
> >> >> > I'm not fitting anything directly.
> >> >> > There is supposed to be a closed form expression, the estimated
> >> >> > coefficient of each polynomial is just the mean of that polynomial
> >> >> > evaluated at the data. The theory and the descriptions sounded
> easy,
> >> >> > but unfortunately it didn't work out.
> >> >> >
> >> >> > I was just hoping to get lucky and that I'm able to skip the small
> >> >> > print.
> >> >> > http://onlinelibrary.wiley.com/doi/10.1002/wics.97/abstract
> >> >> > got me started and it has the fourier case that works.
> >> >> >
> >> >> > There are lots of older papers that I only skimmed, but I should be
> >> >> > able to work out the Hermite case before going back to the general
> >> >> > case again with arbitrary orthogonal polynomial bases. (Initially I
> >> >> > wanted to ignore the Hermite bases because gaussian_kde works well
> in
> >> >> > that case.)
> >> >>
> >> >> Just another graph before stopping with this
> >> >>
> >> >> chebyt work if I cheat (rescale at the end
> >> >>
> >> >> f_hat = (f_hat - f_hat.min())
> >> >> fint2 = integrate.trapz(f_hat, grid)
> >> >> f_hat /= fint2
> >> >>
> >> >> graph is with chebyt with 30 polynomials after shifting and scaling,
> >> >> 20 polynomials looks also good.
> >> >>
> >> >> (hunting for the missing scaling term is left for some other day)
> >> >>
> >> >> In any case, there's a recipe for non-parametric density estimation
> >> >> with compact support.
> >> >>
> >> >
> >> > Ah, now I see what is going on -- monte carlo integration to get the
> >> > expansion of the pdf in terms of orthogonal polynomials. So yes, I
> think
> >> > Lagrange polynomials are probably the ones to use unless you use the
> >> > weight
> >> > in the integral. Note that 1.6 also has the Hermite and Laguerre
> >> > polynomials. But it seems that for these things it would also be
> >> > desirable
> >> > to have the normalization constants.
> >>
> >
> > Heh, I meant Legendre.
> >
> >>
> >> It's also intended to be used as a density estimator for real data,
> >> but the idea is the same.
> >>
> >> My main problem seems to be that I haven't figured out the
> >> normalization (constants) that I'm supposed to use.
> >
> > The normalization for Legendre functions over an interval of length L
> would
> > be (2/L)*2/(2*i + 1)*(1/n), where i is the degree of the polynomial, and
> n
> > is the number of samples.
> >
> >>
> >> Given the Wikipedia page that Andy pointed out, I added an additional
> >> weighting term. I still need to shift and rescale, but the graphs look
> >> nice. (the last one, promised) chebyt with 30 polynomials on sample
> >> with 10000 observations.
> >>
> >> (I don't know yet if I want to use the new polynomials in numpy even
> >> thought they are much nicer. I just gave up trying to get statsmodels
> >> compatible with numpy 1.3 because I'm using the polynomials introduced
> >> in 1.4)
> >>
> >> (I will also look at Andy's pointer to Gram–Schmidt, because I
> >> realized that for nonlinear trend estimation I want orthogonal for
> >> discretely evaluated points instead of for the integral.)
> >>
> >
> > QR is Gram-Schmidt. You can use any polynomial basis for the columns.
> There
> > is also a method due to Moler to compute the polynomials using the fact
> that
> > they satisfy a three term recursion, but QR is simpler.
> >
> > Anne suggested some time ago that I should include the Gauss points and
> > weights in the polynomial classes and I've come to the conclusion she was
> > right. Looks like I should include the normalization factors also.
>
> Since I'm going through this backwards, program first, figure out what
> I'm doing, I always appreciate any of these premade helper functions.
>
> A mixture between trial and error and reading Wikipedia. This is the
> reweighted chebychev T polynomial, orthonormal with respect to uniform
> integration,
>
> for i,p in enumerate(polys[:5]):
>    for j,p2 in enumerate(polys[:5]):
>        print i,j,integrate.quad(lambda x: p(x)*p2(x), -1,1)[0]
>

> And now the rescaling essentially doesn't have an effect anymore
> f_hat.min() 0.00393264959543
> fint2 1.00082015654   integral of estimated pdf
>
>
> class ChtPoly(object):
>
>    def __init__(self, order):
>        self.order = order
>        from scipy.special import chebyt
>        self.poly = chebyt(order)
>
>    def __call__(self, x):
>        if self.order == 0:
>            return np.ones_like(x) / (1-x**2)**(1/4.) / np.sqrt(np.pi)
>        else:
>            return self.poly(x) / (1-x**2)**(1/4.) / np.sqrt(np.pi) *
> np.sqrt(2)
>
> If I can get the same for some of the other polynomials, I would be happy.
>
>
The virtue of the Legendre polynomials here is that you don't need the
weight to make them orthogonal. For the Chebyshev I'd be tempted to have a
separate weight function w(x) = 1/(1 - x**2)**.5 and do (1/n)\sum_i
T_n(x_i)*w(x_i), giving a result in a normal polynomial series. The
additional normalization factor due to the interval would then be 2/L in
addition to the terms you already have. The singularity in the weight could
be a problem so Legendre polynomials might be a safer choice.

To simplify this a bit, in the Legendre case you can use legvander(x,
deg=n).sum(0) and scale the results with the factors in the preceding post
to get the coefficients. The legvander function is a rather small one and
you could pull it out and modify it to do the sum on the fly if space is a
problem.

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