Computes a median polish of a two-way table.
TABLE NROW by NCOL matrix containing the table. (Input)
MAXIT Maximum
number of polishing iterations to be performed. (Input)
An
iteration is counted each time the rows or the columns are polished. The
iterations begin by polishing the rows.
PTABLE (NROW + 1) by (NCOL + 1) matrix containing the cell residuals from the fitted table and, in the last row and column, the marginal residuals. (Output)
NROW Number of
rows in the table. (Input)
Default: NROW
= size (TABLE,1).
NCOL Number of
columns in the table. (Input)
Default: NCOL
= size (TABLE,2).
LDTABL Leading
dimension of TABLE
exactly as specified in the dimension statement in the calling
program. (Input)
Default: LDTABL
= size (TABLE,1).
LDPTAB Leading
dimension of PTABLE exactly as specified in
the dimension statement in the calling program. (Input)
Default:
LDPTAB
= size (PTABLE,1).
ITER Number of iterations actually performed. (Output)
Generic: CALL MEDPL (TABLE, MAXIT, PTABLE [, ])
Specific: The specific interface names are S_MEDPL and D_MEDPL.
Single: CALL MEDPL (NROW, NCOL, TABLE, LDTABL, MAXIT, PTABLE, LDPTAB, ITER)
Double: The double precision name is DMEDPL.
The routine MEDPL
performs a median polish on a two-way table. It first copies TABLE
into PTABLE and fills the last row and
last column of PTABLE
with zeroes. It then computes the
row-wise medians, adds these to the
values in the last column and corresponding row, and subtracts them from the
other entries in the corresponding row. Similar computations are then performed
for all NCOL + 1 columns. The whole
procedure is then repeated (using NROW
+ 1 rows) until convergence is achieved (until no changes occur), or until MAXIT
iterations are performed. Convergence is known to have occurred if ITER
is less than MAXIT.
As Emerson and Hoaglin (1983) discuss, it is not necessarily desirable to continue until convergence. If MAXIT is set to twice the maximum of the number of rows and columns plus five, it is likely that convergence will occur.
As Emerson and Hoaglin point out, median polish starting
with rows can lead to a different fit from that obtained by starting with
columns. Although MEDPL
does not make provision for choosing which dimension to start with, TABLE
can be transposed by use of routine
TRNRR
(IMSL MATH/LIBRARY). Use of the transposed table in MEDPL
would result in the iterations beginning with the columns of the original table.
Further descriptions of median polish, which was first proposed by John Tukey,
and examples of its use can be found in Tukey (1977, Chapter 11) and in Velleman
and Hoaglin (1981, Chapter 8).
Workspace may be explicitly
provided, if desired, by use of M2DPL/DM2DPL.
The
reference is:
CALL M2DPL (NROW, NCOL, TABLE, LDTABL, MAXIT, PTABLE, LDPTAB, ITER, WK)
The additional argument is:
WK Work vector of length max(NROW, NCOL).
This example is taken from Emerson and Hoaglin (1983, page 168). It involves data on infant mortality in the United States, classified by fathers education and by region of the country. In order to show the difference between making only one polishing pass over the rows and polishing until convergence, on the first invocation MAXIT is set to one. On a second call, it is set large enough to have reasonable assurance of execution until convergence. In the first case, the last row and column of PTABLE are printed. The values in these are the medians before any polishing. These values approach zero as the polishing continues.
USE MEDPL_INT
USE UMACH_INT
USE WRRRL_INT
IMPLICIT NONE
INTEGER NCOL, NROW
PARAMETER (NCOL=5, NROW=4)
!
INTEGER ITER, LDPTAB, LDTABL, MAXIT, NOUT
REAL PTABLE(NROW+1,NCOL+1), TABLE(NROW,NCOL)
CHARACTER CLABEL(1)*5, RLABEL(1)*5
!
DATA CLABEL/'NONE'/
DATA RLABEL/'NONE'/
DATA TABLE/25.3, 32.1, 38.8, 25.4, 25.3, 29.0, 31.0, 21.1, 18.2, &
18.8, 19.3, 20.3, 18.3, 24.3, 15.7, 24.0, 16.3, 19.0, 16.8, &
17.5/
!
CALL UMACH (2, NOUT)
MAXIT = 1
LDTABL = 4
LDPTAB = 5
CALL MEDPL (TABLE, MAXIT, PTABLE, ITER=ITER)
CALL WRRRL ('Fitted table after one iteration over the rows', &
PTABLE, CLABEL, RLABEL, FMT='(W10.4)')
MAXIT = 15
CALL MEDPL (TABLE, MAXIT, PTABLE, ITER=ITER)
CALL WRRRL ('%/Fitted table and marginal residuals', PTABLE, &
CLABEL, RLABEL, FMT='(W10.4)')
WRITE (NOUT,99999) ITER
99999 FORMAT (/, ' Iterations taken: ', I2)
END
Fitted table after one
iteration over the
rows
7.0
7.0
-0.1
0.0
-2.0
18.3
7.8
4.7
-5.5
0.0
-5.3
24.3
19.5
11.7
0.0
-3.6
-2.5
19.3
4.3
0.0
-0.8
2.9
-3.6
21.1
0.0
0.0
0.0
0.0
0.0
0.0
Fitted table and marginal
residuals
-1.55
0.00
0.00
-1.15
0.60
-1.45
1.55
0.00
-3.10
1.15
-0.40
2.25
10.85
4.60
0.00
-4.85
0.00
-0.35
-3.25
-6.00
0.30
2.75
0.00
0.35
8.10
6.55
-0.55
0.70
-3.05 20.20
Iterations taken: 15
PHONE: 713.784.3131 FAX:713.781.9260 |