Chapter 2: Regression

OPOLY

Generates orthogonal polynomials with respect to x-values and specified weights.

Required Arguments

X — Vector of length N containing the x-values.   (Input)

NDEG — Degree of highest degree orthogonal polynomial to be generated.   (Input)

SMULTC — Multiplicative constant used to compute a scaled version of x on the interval 2 to 2, inclusive.   (Output)

SADDC — Additive constant used to compute a scaled version of x on the
interval 2 to 2, inclusive.   (Output)

SX — Vector of length N containing the scaled version of x on the interval 2 to 2, inclusive, computed as follows: SX(i) = SMULTC * X(i) + SADDC where i = 1, 2, …, N.   (Output)
If X is not needed, SX and X can occupy the same storage locations.

A — Vector of length NDEG containing constants used to generate orthogonal polynomials.   (Output)

B — Vector of length NDEG containing constants used to generate orthogonal polynomials.   (Output)

POLY — Matrix, N by NDEG, containing the orthogonal polynomials evaluated at SX(i) for
i = 1, 2, …, N.   (Output)

Optional Arguments

N — Number of x-values.   (Input)
Default: N = size (POLY,1).

IWT — Weighting option.   (Input)
IWT = 0 means that all weights are 1.0. For IWT = 1, WT contains the weights.
Default: IWT = 0.

WT — Vector of length N containing the weights.   (Input, if IWT = 1)
If IWT = 0, WT is not referenced and can be a vector of length one.

LDPOLY — Leading dimension of POLY exactly as specified in the dimension statement in the calling program.   (Input)
Default: LDPOLY = size (POLY,1).

FORTRAN 90 Interface

Generic:                              CALL OPOLY (X, NDEG, SMULTC, SADDC, SX, A, B, POLY [,…])

Specific:                             The specific interface names are S_OPOLY and D_OPOLY.

FORTRAN 77 Interface

Single:            CALL OPOLY (N, X, IWT, WT, NDEG, SMULTC, SADDC, SX, A, B, POLY, LDPOLY)

Double:                              The double precision name is DOPOLY.

Description

Routine OPOLY generates orthogonal polynomials over a set of xi’s and with respect to weights wi. The routine OPOLY is based on the algorithm of Forsythe (1957). (See also Kennedy and Gentle 1980.) A modification to Forsythe’s algorithm is made for the inclusion of weights (Kelly 1967, page 68).

Let xi be a value of the independent variable. The xi’s are scaled to the interval [2, 2] for computational accuracy. The scaled version of the independent variable is computed by the formula zi= mxi+ c. The multiplicative scaling constant m (stored in SMULTC) is

The additive constant c (stored in SADDC) is

Orthogonal polynomials are generated using the three-term recurrence relationship

beginning with the initial polynomials

The aj’s and bj’s (stored in A and B) are computed to make the pj(z)’s orthogonal, with respect to the the set of weights wi, and over the set zi.

Comments

1.         Informational error

Type Code

3         8                  N must be greater than NDEG in order for higher order polynomials to be nonzero. Columns N + 1 through NDEG of POLY are set to zero.

2.         The orthogonal polynomials evaluated at each scaled X value are computed from A and B as follows:
POLY(I, 1) = SX(I) A(1)
POLY(I, 2) = (SX(I) A(2)) * POLY(I, 1) B(2)
POLY(I, J) = (SX(I) A(J)) * POLY(I, J 1) B(J) * POLY(I, J 2) for J = 3 through NDEG.

3.         If NDEG is greater than 10, the accuracy of the results may be questionable.

Example

First-degree and second-degree orthogonal polynomials are generated using equally spaced x values 1, 2, …, 12. (Equally spaced x values are not required by OPOLY.)

 

      USE OPOLY_INT

      USE UMACH_INT

      USE WRRRN_INT

 

      IMPLICIT   NONE

      INTEGER    LDPOLY, N, NDEG

      PARAMETER  (N=12, NDEG=2, LDPOLY=N)

!

      INTEGER    NOUT

      REAL       A(NDEG), B(NDEG), POLY(LDPOLY,NDEG), SADDC, SMULTC, &

                 SX(N), WT(1), X(N)

!

      DATA X/1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, &

           12.0/

!

      CALL OPOLY (X, NDEG, SMULTC, SADDC, SX, A, B, POLY)

      CALL UMACH (2, NOUT)

      WRITE (NOUT,99999) SMULTC, SADDC

99999 FORMAT (' SMULTC = ', F7.3, '  SADDC = ', F7.3)

      CALL WRRRN ('A', A, 1, NDEG, 1)

      CALL WRRRN ('B', B, 1, NDEG, 1)

      CALL WRRRN ('POLY', POLY)

      END

Output

 

SMULTC =   0.364  SADDC =  -2.364

           A
         1           2
-5.960E-08  -1.009E-07

      B
    1       2
0.000   1.576

       POLY
         1       2
 1  -2.000   2.424
 2  -1.636   1.102
 3  -1.273   0.044
 4  -0.909  -0.749
 5  -0.545  -1.278
 6  -0.182  -1.543
 7   0.182  -1.543
 8   0.545  -1.278
 9   0.909  -0.749
10   1.273   0.044
11   1.636   1.102
12   2.000   2.424



http://www.vni.com/
PHONE: 713.784.3131
FAX:713.781.9260