[SciPy-User] orthonormal polynomials (cont.)

[josef.pktd at gmail.com]

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)

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.

Josef

>
> <snip>
>
> Chuck
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
>