[SciPy-User] orthonormal polynomials (cont.)

[josef.pktd at gmail.com]

On Thu, Jun 16, 2011 at 1:42 PM, Charles R Harris
<charlesr.harris at gmail.com> wrote:
>
>
> 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.

For me this is still scipy.special, which is not my area, and not yet stats.

This is just to construct a polynomial basis, that then can be used
for density estimation or testing, similar to the last time I did
this. (I have not gotten this far yet with the new version, I'm still
writing tests for emerson.)

we can get a density estimate with

pdf_estimated(x) = w(x)  sum_{i}  ahat_{i} * poly_{i} (x)

with ahat{i} = sum_{k} poly_{i} (xdata_{k})     that's the same as
last time, and where the actual data comes in.

w(x) is a base density, and the rest is a polynomial expansion around
it. In my previous examples, I assumed w(x) is 1 (Lebesque measure)

There are some variation on how the expansion of the density is done,
but all start with a polynomial basis for the expansion.

I hope to have some examples soon.

thanks,

Josef

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