[SciPy-User] orthonormal polynomials (cont.)

[Charles R Harris]

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 ;)

<snip>

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