[SciPy-User] orthonormal polynomials (cont.)

[Charles R Harris]

On Thu, Jun 16, 2011 at 10:38 AM, <josef.pktd at gmail.com> wrote:

> On Thu, Jun 16, 2011 at 12:00 PM, Charles R Harris
> <charlesr.harris at gmail.com> wrote:
> >
> >
> > On Thu, Jun 16, 2011 at 9:29 AM, <josef.pktd at gmail.com> wrote:
> >>
> >> .......... <still working>
> >>
> >> example with Emerson recursion
> >>
> >> laguerre for weight function gamma pdf with shape=0.5 instead of shape=1
> >>
> >> >>> mnc = moments_gamma(np.arange(21), 0.5, 1)
> >> >>> myp = orthonormalpoly_moments(mnc, 10, scale=1)
> >> >>> innerprod = dop.inner_cont(myp, 0., 100, stats.gamma(0.5).pdf)[0]
> >> >>> np.max(np.abs(innerprod - np.eye(innerprod.shape[0])))
> >> 5.1360515840315202e-10
> >> >>> for pi in myp: print pi.coeffs
> >> ...
> >> [ 1.]
> >> [ 1.41421356 -0.70710678]
> >> [ 0.81649658 -2.44948974  0.61237244]
> >> [ 0.2981424  -2.23606798  3.35410197 -0.55901699]
> >> [ 0.07968191 -1.1155467   4.18330013 -4.18330013  0.52291252]
> >> [ 0.01679842 -0.37796447  2.64575131 -6.61437828  4.96078371
> -0.49607837]
> >> [  2.92422976e-03  -9.64995819e-02   1.08562030e+00  -5.06622805e+00
> >>   9.49917760e+00  -5.69950656e+00   4.74958880e-01]
> >> [  4.33516662e-04  -1.97250081e-02   3.25462634e-01  -2.44096975e+00
> >>   8.54339413e+00  -1.28150912e+01   6.40754560e+00  -4.57681829e-01]
> >> [  5.59667604e-05  -3.35800562e-03   7.63946279e-02  -8.40340907e-01
> >>   4.72691760e+00  -1.32353693e+01   1.65442116e+01  -7.09037640e+00
> >>   4.43148525e-01]
> >> [  6.39881348e-06  -4.89509231e-04   1.46852769e-02  -2.22726700e-01
> >>   1.83749528e+00  -8.26872874e+00   1.92937004e+01  -2.06718219e+01
> >>   7.75193320e+00  -4.30662955e-01]
> >>
> >> I haven't had time to figure out QR yet.
> >>
> >
> > One way to think of QR in this context is to think of the columns of the
> > Vandermonde matrix V as the basis functions. Then
> >
> > QR = V
> > Q = V*R^-1
> >
> > Since R and its inverse are upper triangular, the orthonormal columns of
> Q
> > are expressed as linear combinations of the basis functions in the
> columns
> > of V of degree <= the column index. For general numerical reasons I would
> > use a Chebyshev basis rather than powers of x.
> >
> > I can't find a reference on the Emerson recursion. I'm guessing that it
> is
> > for a power series basis and generates new columns on the fly as c_i =
> > x*c_{i -1} so that the new column is orthogonal to all the columns c_j, j
> <
> > i - 2. Anyway, that's what I would do if I wanted better numerical
> > conditioning ;)
>
> I pretty much implemented this after I discovered that they use
> zero-based indexing and it didn't look too scary
>
>    J. C. W Rayner, O. Thas, and B. De Boeck, “A GENERALIZED EMERSON
> RECURRENCE
>    RELATION,” Australian & New Zealand Journal of Statistics 50, no. 3
>    (September 1, 2008): 235-240.
>
> http://onlinelibrary.wiley.com/doi/10.1111/j.1467-842X.2008.00514.x/abstract
>
> Before that, I tried your QR example for a while but I had two problems,
>
> * the resulting polynomials are not orthogonal if I integrate,  int
> poly_i(x) * poly_j(x) dx
> * I need orthogonality with respect to a weight function:  int
> poly_i(x) * poly_j(x) * w(x) dx == (i==j).astype(int)
>
>
I think the question is how you perform the integration. The QR does it
numerically with the sample points passed into the *vander functions and I
used uniform spacing for uniform measure. Weighting the rows with the sqrt
of the weight function will produce polynomials orthogonal for that weight.
The whole thing can be vastly improved by using selected sample points if
the weight function is an actual function that can be evaluated at arbitrary
points. Send me an example and I will work it out for you.

The first I may not need anymore. emerson works for continuous functions.
> The second I would like to figure out when I move to discrete
> distribution, where I have sum instead of integral. (But after I
> finish with the continuous distributions).
> sum_{x in X}  poly_i(x) * poly_j(x) * w(x) dx == (i==j).astype(int)
> Is there a way to get weights into QR?
>
> The Emerson recursion that I have only works with power series, so I
> still don't know how to do it with any other basis functions.
>
>
If it is what I think it is it shouldn't be difficult. I can't get to the
reference you linked.

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