MEDPL
Computes a median polish of a two-way table.
Required Arguments
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)
Optional Arguments
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)
FORTRAN 90 Interface
Generic: CALL MEDPL (TABLE, MAXIT, PTABLE [, …])
Specific: The specific interface names are S_MEDPL and D_MEDPL.
FORTRAN 77 Interface
Single: CALL MEDPL (NROW, NCOL, TABLE, LDTABL, MAXIT, PTABLE, LDPTAB, ITER)
Double: The double precision name is DMEDPL.
Description
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).
Comments
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).
Example
This example is taken from Emerson and Hoaglin (1983, page 168). It involves data on infant mortality in the United States, classified by father’s 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
Output
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
