[SciPy-User] orthonormal polynomials (cont.)

[Charles R Harris]

On Thu, Jun 16, 2011 at 11:32 AM, Charles R Harris <
charlesr.harris at gmail.com> wrote:

>
>
> 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.
>
>
OK, I found it, no surprises, it's the standard three term recurrence with
expectations replacing integrals. Are you using sampled data? I thought you
wanted polynomials for a specified weight over an interval.

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