(FILECREATED "16-Feb-86 13:03:43" {QV}<IDL>SOURCES>COMPRESS.;16 41852
changes to: (FNS COVAR PAIRN)
previous date: "27-Nov-85 17:40:38" {QV}<IDL>SOURCES>COMPRESS.;15)
(* Copyright (c) 1983, 1984, 1985, 1986 by Xerox Corporation. All rights reserved.)
(PRETTYCOMPRINT COMPRESSCOMS)
(RPAQQ COMPRESSCOMS ((* Data matrix compression routines)
(FNS COUNTS COVAR GROUP GROUP.HIDDENCODES GROUP.HIDDENLEVLAB GROUP.LAB
GROUP.LINLOC GROUP.SHAPE MOMENTS MOMENTS.FEW MOMENTS.MANY PAIRN POOL)))
(* Data matrix compression routines)
(DEFINEQ
(COUNTS
[ULAMBDA ((A (EXPECTS ARRAY))
(RETURNS ARITH))
(* bas: "15-FEB-83 17:48")
(* User entry to counts -
i.e. summation without NIL producing NIL)
(if (IEQP 0 (NLOGICALELTS (fetch SHAPE of A)))
then (if (VSCALARP A)
then (COPYAELT A (VSCALARPTR A))
else 0)
else (bind V (TOT ← 0)
(GSBA ←(SETUP A (QUOTE ROWMAJOR)))
declare (V SCALAR)
(TOT ARITH)
(GSBA GENSTATEBLOCK)
until (fetch DONE of GSBA) when (SETQ V (GETAELT A (NEXT GSBA)))
do (SETQ TOT (PLUS TOT V)) finally (RETURN TOT)))])
(COVAR
[ULAMBDA ((A (EXPECTS MATRIX))
(WT (EXPECTS VECTOR (ONEOF NIL VECTOR)) (* Weight vector))
(RETURNS MATRIX))
(* edited: "16-Feb-86 12:38" posted: "27-MAY-77 23:27"
) (* This is the covar product compression routine.
It always returns the covar matrix, which will have
NIL's for every pair of disjoint variables.)
(COND
((AND WT (NEQ (GETRELT (fetch SHAPE of WT)
1)
(GETRELT (fetch SHAPE of A)
1)))
(UERROR "Wrong length weighting vector")))
(DPROG ((COVAR NIL SIMARRAY (* Output matrix))
(NVARS (GETRELT (fetch SHAPE of A)
2) IJK (* Number of variables))
(MFLAG NIL BOOL)
(N 0.0 FLOATING (SATISFIES (LEQ 0.0 N)) (* Sum of the weights))
(NROW NIL ROWFLOAT (* If needed, pairwise N's))
(MJ NIL ROWFLOAT (* Pairwise means))
(MK NIL ROWFLOAT)
(SJ NIL ROWFLOAT (* Pairwise sums of squares))
(SK NIL ROWFLOAT)
(GSBW (AND WT (SETUP WT (QUOTE ROWMAJOR))) (ONEOF NIL GENSTATEBLOCK))
(GSBA (SETUP A (QUOTE ROWMAJOR)) GENSTATEBLOCK)
THEN (OBS (create ROWFLOAT
NELTS ← NVARS) ROWSCALAR)
(NR (IQUOTIENT (ITIMES NVARS (ADD1 NVARS))
2) IJK)
THEN (FIRSTMEAN (ADD1 NR) IJK (* Index in COVROW of the first mean))
(LASTMEAN (IPLUS NR NVARS) IJK (* Index in COVROW of the last mean))
THEN (COVROW (create ROWFLOAT
NELTS ←(ADD1 LASTMEAN)
INIT ← 0.0) ROWFLOAT (* Covar accumulation. The ADD1 allows for the -1/N)))
(* COVROW is treated as a row for most of this routine, but it then becomes the elementblock of the output array.
Thus, we make certain assumptions about the correspondence between the layout of element blocks and the addressing
on rows. There is a further trick: the mean estimates are accumulated in the last row of COVROW, where they
eventually will end up. See the comment on MEANJ below)
[bind (WVAL ← 1.0) until (fetch DONE of GSBA)
do
(COND
((NULL WT))
((AND (SETQ WVAL (GETAELT WT (NEXT GSBW)))
(FGREATERP (SETQ WVAL (FLOAT WVAL))
0.0)))
(T (SKIP GSBA NVARS) (* Skip this whole subject)
(GO $$ITERATE))) (* Read observations for each var)
[for J to NVARS unless (SETRELT OBS J (GETAELT A (NEXT GSBA)))
do (COND
((NOT MFLAG)
(SETQ MFLAG T) (* The very first missing observation)
(SETQ NROW (create ROWFLOAT
NELTS ← NR
INIT ← N))
(SETQ MK (create ROWFLOAT
NELTS ← NR))
(SETQ MJ (create ROWFLOAT
NELTS ← NR))
(SETQ SJ (create ROWFLOAT
NELTS ← NR))
(SETQ SK (create ROWFLOAT
NELTS ← NR))
(for I (IK ← 0) to NVARS as IM from FIRSTMEAN as DI from 1
by I
declare (IK IJK)
(IM IJK)
(DI IJK)
do (for K to I as KM from FIRSTMEAN as DK from 1
by K
declare (KM IJK)
(DK IJK)
do (add IK 1)
(SETRELT MJ IK (GETRELT COVROW IM))
(SETRELT MK IK (GETRELT COVROW KM))
(SETRELT SJ IK (GETRELT COVROW DI))
(SETRELT SK IK (GETRELT COVROW DK]
(bind DJ DK OBSJ P NSC (JK ← 0)
declare (DJ FLOATING)
(DK (ONEOF NIL FLOATING))
(JK IJK)
(OBSJ (ONEOF NIL FLOATING))
(P FLOATING (* Scratch value))
(NSC FLOATING (* N Scaling factor = (N+W) *N/W))
(WVAL FLOATING)
for J to NVARS as MEANJ from FIRSTMEAN declare (MEANJ IJK)
first [COND
((NOT MFLAG) (* If no missing data set the new N and the N
adjustment factor outside the loops)
(SETQ NSC (FQUOTIENT (FTIMES N (SETQ N (FPLUS N WVAL)))
WVAL]
do
(COND
((NULL (SETQ OBSJ (GETRELT OBS J)))
(add JK J) (* Skip that row of the output)
)
(T (DPROGN ((OBSJ FLOATING))
(COND
((NOT MFLAG)
(* If no missing data the observation is converted into a mean adjustment and stored back into the OBS row in that
form. This is done so that the K loop can pick up the value without having to rescale it, as the N on which the
scaling depends is a constant throughout.)
[SETRELT OBS J (SETQ DJ (FTIMES (FQUOTIENT WVAL N)
(FDIFFERENCE OBSJ
(GETRELT COVROW
MEANJ]
(fadd (GETRELT COVROW MEANJ)
DJ)))
(for K to J
do (add JK 1) (* Move to the next cell of the output)
(COND
((NULL (SETQ DK (GETRELT OBS K)))
(GO $$ITERATE)))
(DPROGN ((DK FLOATING))
[COND
(MFLAG [SETRELT NROW JK (SETQ N (FPLUS WVAL
(SETQ NSC
(GETRELT
NROW JK]
(* Update the pairwise N)
(SETQ NSC (FQUOTIENT (FTIMES N NSC)
WVAL))
(* Set the N adjustment factor for this N)
(SETQ P (FQUOTIENT WVAL N))
(* Weighting factor for this subject)
(COND
((EQ K J)
[SETQ DK (SETQ DJ
(FTIMES P (FDIFFERENCE OBSJ
(GETRELT COVROW
MEANJ]
(fadd (GETRELT COVROW MEANJ)
DJ)
(* Don't use the diagonal of MJ MK SJ or SK -
the results are elsewhere)
)
(T [fadd (GETRELT MJ JK)
(SETQ DJ (FTIMES P
(FDIFFERENCE
OBSJ
(GETRELT MJ JK]
[fadd (GETRELT MK JK)
(SETQ DK (FTIMES P
(FDIFFERENCE
DK
(GETRELT MK JK]
(* Update the pairwise means)
(fadd (GETRELT SJ JK)
(FTIMES DJ DJ NSC))
(fadd (GETRELT SK JK)
(FTIMES DK DK NSC))
(* Update the pairwise sums of squares)
]
(fadd (GETRELT COVROW JK)
(FTIMES DJ DK NSC))
(* No matter which path you are on, update the
crossproducts)
)))]
[SETQ COVAR (create SIMARRAY
FORMAT ←(QUOTE SYMMETRIC)
ELEMENTBLOCK ← COVROW
SHAPE ←(create ROWINT
NELTS ← 2
INIT ←(ADD1 NVARS]
(COND
(LABPROPFLAG [SETTITLE COVAR (MAKETITLE (QUOTE COVARIATIONS)
A
(COND
(WT (LIST ", Weighted by " WT]
(LAB.COPYDIM A COVAR 2 1)
(SETLEVLAB COVAR 1 (ADD1 NVARS)
(QUOTE Constant))
(LAB.COPYDIM COVAR COVAR 1 2)))
[COND
(MFLAG (SETARRAYPROP COVAR (QUOTE PAIRN)
T)
(bind TEMP declare (TEMP FLOATING) for K to NR
when (AND (FGREATERP N (SETQ TEMP (GETRELT NROW K)))
(FGREATERP TEMP 0.0))
do (SETQ N TEMP)) (* Setting N to the smallest pairwise N)
(COND
((FGTP N 0.0)
(bind TEMP PD (RC ← 0)
(CEB ←(fetch ELEMENTBLOCK of COVAR))
declare (TEMP FLOATING)
(PD FLOATING (* Row diagonal elt))
(RC IJK)
(CEB ROWSCALAR (* Loosens COVROW in a formally illegal way)
)
for R to NVARS as D from 1 by R declare (D IJK)
do [SETQ PD (SETRELT COVROW D (FTIMES (GETRELT COVROW D)
(FQUOTIENT N
(GETRELT NROW D]
(* Adjusting the sums of squares to the minimum
pairwise N)
[for C to (SUB1 R) as D from 1 by C declare (D IJK)
do (add RC 1)
(SETRELT CEB RC
(COND
((FGTP (SETQ TEMP (GETRELT NROW RC))
0.0)
(* Correcting the crossproduct to be the projection of
the scaled sums of squares by the correlation
coefficient computed on all pairwise present data)
(FTIMES (GETRELT COVROW RC)
(SQRT (FQUOTIENT
(FTIMES PD
(GETRELT COVROW
D))
(FTIMES (GETRELT SJ RC)
(GETRELT SK RC]
(add RC 1) (* Skip the diagonal entry--it's already been done)
]
(COND
((FGTP N 0.0)
(SETRELT COVROW (ADD1 LASTMEAN)
(FQUOTIENT -1.0 N)))
(T (DPROG ((CEB (fetch ELEMENTBLOCK of COVAR) ROWSCALAR))
(ROWBLT CEB NIL)) (* Only happens for completely missing or empty
matrix)
))
(RETURN COVAR))])
(GROUP
[ULAMBDA ((ATTRIBS (ONEOF VECTOR MATRIX))
(VALUES (ONEOF NIL SCALAR ARRAY))
(DIM ANY (* The nested dimensions of VALUES, 1 if NIL)))
(* jop: "27-Nov-85 17:36")
(if (NULL VALUES)
then (SETQ VALUES 1))
(DPROG ((ASHP (fetch SHAPE of ATTRIBS) ROWINT)
THEN (NOBS (GETRELT ASHP 1) IJK))
(if (type? ARRAY VALUES)
then (SETQ DIM (MAKE1DIMSPEC VALUES (OR DIM 1)))
(OR (IEQP (GETRELT (fetch SHAPE of VALUES)
DIM)
NOBS)
(UERROR "Attributes-Values length mismatch"))
else (SETQ VALUES (VFROMR (create ROWINT
NELTS ← NOBS
INIT ← VALUES)
(MAKETITLE VALUES)))
[coerce DIM (ONEOF NIL (INTEGER (SATISFIES (EQ DIM 1]
(* Just to check the user)
(SETQ DIM 1))
(RETURN
(DPROG ((GROUP NIL SIMARRAY)
(VSHP (fetch SHAPE of VALUES) ROWINT)
(NDV (fetch NDIMS of VALUES) INTEGER)
(NFAC (if (EQ (fetch NELTS of ASHP)
1)
then 1
else (GETRELT ASHP 2)) INTEGER)
THEN (CBROW (create ROWPTR
NELTS ← NFAC) ROWPTR (* Row of codebooks))
(GDIM (IPLUS NFAC DIM) INTEGER))
[if (ZEROP NFAC)
then (* If there are no factors, then simply return a copy
of VALUES, which will be be appropriately labeled and
titled below.)
(SETQ GROUP (PRESERVE VALUES))
else
(DPROG ((ACC NIL ROWINT (* Number of items per cell, then initial storage
locations))
(DPOS (create ROWINT
NELTS ← NOBS
INIT ← 0) ROWINT
(* Maps items into their classifier cells))
(MAXCNT 0 INTEGER (* Largest cell N))
(GRSHP (create ROWINT
NELTS ←(IPLUS NFAC NDV)) ROWINT))
(GROUP.SHAPE ATTRIBS GRSHP CBROW)
(for J from 1 to NDV do (SETRELT GRSHP (IPLUS NFAC J)
(GETRELT VSHP J)))
(* Fill in shape for trailing dimensions from VALUES)
(for I (NCELLS ← 1) to NFAC
declare (I INTEGER)
(NCELLS IJK)
do [if (IEQP (SETQ NCELLS (ITIMES NCELLS (GETRELT GRSHP I)))
0)
then (* Pack together some garbage for the silly user)
(RETURN (ALLOC.SARRAY GRSHP (QUOTE FULL]
finally (SETQ ACC (create ROWINT
NELTS ← NCELLS
INIT ← 0)))
(DPROG ((OSET (create ROWINT
NELTS ← NFAC) ROWINT
(* Quasi offsets for GROUP.LINLOC)))
(SETRELT OSET NFAC 1)
[for I (V ← 1) from NFAC to 2 by -1
do (SETQ V (SETRELT OSET (SUB1 I)
(ITIMES V (GETRELT GRSHP I]
(for I to NOBS bind P (GSB ←(SETUP ATTRIBS (QUOTE ROWMAJOR)))
declare (P (ONEOF NIL IJK)
(* Location in the classifier space)
)
(GSB GENSTATEBLOCK)
do (if (SETQ P (GROUP.LINLOC ATTRIBS NFAC GSB OSET CBROW))
then (SETQ MAXCNT (IMAX MAXCNT
(add (GETRELT ACC P)
1)))
(SETRELT DPOS I P))))
(SETRELT GRSHP GDIM MAXCNT)
(SETQ GROUP (create SIMARRAY
SHAPE ← GRSHP
FORMAT ←(QUOTE FULL)
ELEMENTBLOCK ←(create ROWSCALAR
NELTS ←(NLOGICALELTS GRSHP)
RELTTYPE ←(fetch AELTTYPE
of VALUES)
INIT ← NIL)))
[bind GSBV VSLICE GDIMINC (WIDTH ←(GETRELT (fetch (SIMARRAY OFFSETS)
of GROUP)
GDIM))
(GSBG ←(SETUP GROUP (QUOTE ROWMAJOR)))
declare (GSBG GENSTATEBLOCK)
(GDIMINC IJK)
(WIDTH IJK)
(VSLICE ARRAY)
(GSBV GENSTATEBLOCK)
first (SETQ GDIMINC (IDIFFERENCE (GETRELT (fetch (SIMARRAY
OFFSETS)
of GROUP)
(SUB1 GDIM))
WIDTH))
[if (EQ NDV 1)
then (* A vector VALUES would yield scalars on selection,
hence is treated specially.)
(SETQ GSBV (SETUP VALUES (QUOTE ROWMAJOR)))
else
(SETQ VSLICE
(DPROG ((LABPROPFLAG NIL SPECIAL))
[RETURN (FSELECT VALUES
(DPROG ((S (create ROWPTR
NELTS ← NDV
INIT ←(QUOTE
ALL)) ROWPTR)
(RETURNS ROWPTR))
(SETRELT S DIM 1)
(RETURN S))])]
(bind (OFF ←(GETRELT (fetch (SIMARRAY OFFSETS) of GROUP)
NFAC))
for J to (fetch NELTS of ACC) as LOC from 0
by OFF
declare (LOC IJK)
(OFF IJK)
do (SETRELT ACC J LOC))
(* Setting up the initial storing locations in the
output)
for I CELL to NOBS
declare (I IJK)
(CELL IJK)
unless (IEQP (SETQ CELL (GETRELT DPOS I))
0)
do (RESETUP GSBG)
(SKIP GSBG (GETRELT ACC CELL))
(add (GETRELT ACC CELL)
WIDTH)
(if (EQ NDV 1)
then (SETAELT GROUP (NEXT GSBG)
(GETAELT VALUES (NEXT GSBV)))
else (ADJUST.SELECTION VSLICE DIM I)
(SETQ GSBV (SETUP VSLICE (QUOTE ROWMAJOR)
GSBV))
(until (fetch DONE of GSBV)
do [for J to WIDTH
do (SETAELT GROUP (NEXT GSBG)
(GETAELT VSLICE (NEXT GSBV]
(* Have to skip cause DIM might not be 1)
(SKIP GSBG GDIMINC])]
[if LABPROPFLAG
then (GROUP.LAB ATTRIBS GROUP CBROW)
(SETDIMLAB GROUP GDIM (OR (GETDIMLAB VALUES DIM)
(GETDIMLAB ATTRIBS 1)
(QUOTE Within)))
(if (IGREATERP (LABELLEVEL VALUES)
1)
then (for I from (ADD1 NFAC) as J from 1 to NDV
unless (EQ J DIM)
do (LAB.COPYDIM VALUES GROUP J I T)))
(SETTITLE GROUP (MAKETITLE NIL VALUES (LIST " by " ATTRIBS]
[replace KEEPS of GROUP with (NCONC (for I to NFAC collect I)
(for I in (fetch KEEPS
of VALUES)
collect
(* Must correct for new dimensions)
(IPLUS NFAC I]
(RETURN GROUP))))])
(GROUP.HIDDENCODES
[DLAMBDA ((A ARRAY)
(RETURNS (ONEOF NIL CODEBOOK)))
(* rmk: " 4-JUL-79 22:20" posted: " 4-JUL-79 22:20")
(* Finds a codebook hidden behind an integer selection
on the valdim)
[AND (type? SELARRAY A)
(DPROG ((TT (A:TTAB) ROWPTR)
(AB (A:BASEARRAY) SIMARRAY)
THEN (VD (GETVALDIM AB) (ONEOF NIL INTEGER))
(RETURNS (ONEOF NIL CODEBOOK)))
(RETURN (if (AND VD (type? INTEGER (TT$VD)))
then (GETCODES AB TT$VD))))]])
(GROUP.HIDDENLEVLAB
[DLAMBDA ((A ARRAY)
(RETURNS LABEL))
(* bas: " 1-DEC-78 16:24" posted: " 1-DEC-78 16:48")
(* Looks for a single integer selection and returns the
hidden level label from such)
(AND (type? SELARRAY A)
(DPROG ((AB (A:BASEARRAY) SIMARRAY)
(TT (A:TTAB) ROWPTR)
(TTE NIL TTELT)
(FOUND NIL BOOL)
(LAB NIL LABEL)
(RETURNS LABEL))
(for I to TT:NELTS
do (TTE←TT$I)
(SELECTQ (TTELTTYPE TTE)
[INTEGER (if FOUND
then (RETURN)
else (FOUND←T)
(LAB←(GETLEVLAB AB I TTE]
(ARRAY (RETURN))
NIL))
(RETURN LAB)))])
(GROUP.LAB
[DLAMBDA ((A (ONEOF VECTOR MATRIX))
(OUTPUT ARRAY)
(CBROW ROWPTR))
(* jop: " 7-Oct-85 23:21")
(* Adds classification labels and title to the
unlabeled OUTPUT array, based on the labels of A and
the design CBROW)
[bind (SHAPE ←(fetch SHAPE of OUTPUT)) declare (SHAPE ROWINT) for DIM
to (fetch NELTS of CBROW)
do [SETDIMLAB OUTPUT DIM (if (EQ (fetch NDIMS of A)
2)
then (OR (GETLEVLAB A 2 DIM)
(PACK* (QUOTE Factor)
DIM))
else (OR (GROUP.HIDDENLEVLAB A)
(QUOTE Value]
(bind (CB ←(GETRELT CBROW DIM)) for LEV to (GETRELT SHAPE DIM)
declare (CB (ONEOF NIL CODEBOOK)) do (SETLEVLAB OUTPUT DIM LEV (perform
CODEBOOK.LABFROMINDEX CB
LEV]])
(GROUP.LINLOC
[DLAMBDA ((A (ONEOF VECTOR MATRIX))
(NFAC INTEGER)
(GSB GENSTATEBLOCK)
(OFFS ROWINT)
(CBROW ROWPTR)
(RETURNS (ONEOF NIL IJK)))
(* bas: "15-FEB-83 17:35")
(* Produces an linloc for SHP for the classification of A represented by CBROW for the next NFAC observations
generated by GSBA. Returns NIL if one of the classifying observations is wild or missing.)
(for VAR IDX (LOC ← 1) to NFAC
declare (IDX (ONEOF NIL INTEGER))
(LOC IJK)
do (if [SETQ IDX (perform CODEBOOK.INDEX (GETRELT CBROW VAR)
(GETAELT A (NEXT GSB]
then (add LOC (ITIMES (GETRELT OFFS VAR)
(SUB1 IDX)))
else (SKIP GSB (IDIFFERENCE NFAC VAR))
(RETURN))
finally (RETURN LOC))])
(GROUP.SHAPE
[DLAMBDA ((A (ONEOF VECTOR MATRIX))
(SHAPE ROWINT (* The shape whose factor part is to be filled in))
(CBROW ROWPTR (* The codebook row to be filled in)))
(* rmk: " 6-JUL-79 12:04" posted: "13-JUL-77 21:52")
(* Computes the shape of an array appropriate for the classification design represented by A.
-
On return, CBROW CAR is the codebook. For variables without a codebook, one is constructed containing all the
distinct values)
(DPROG ((NVARS (CBROW:NELTS) INTEGER)
(STARTPTR NIL INTEGER)
(NOCBCOUNT 0 INTEGER))
(* Determines the non-codebook columns and builds in SHAPE a threaded list of these columns.
If for example the columns 2 3 5 and 6 of an array with six columns are non-codebook, then SHAPE will be {FOO 3 5
FIE 6 2}. startptr will be 2- the first non-codebook column. The position in SHAPE indicates the corresponding
position in CBROW where the codebook is constructed. In the example, FOO and FIE stand for the length of the
codebooks for codebook-present levels. They are inserted on the first past through.)
(DPROG ((CODES NIL (ONEOF NIL CODEBOOK)))
(SELECTQ (GETVALDIM A)
[2 (for I PTR to NVARS do (if CODES←(GETCODES A I)
then ((SHAPE$I)←(perform CODEBOOK.LENGTH CODES))
(CBROW$I←CODES)
else (if PTR
then ((SHAPE$PTR)←I)
else STARTPTR←I)
(PTR←I)
(add NOCBCOUNT 1))
finally (if PTR
then ((SHAPE$PTR)←STARTPTR]
(1 (NOCBCOUNT←NVARS)
((SHAPE$NVARS)←STARTPTR←1)
(for I to NVARS-1 do ((SHAPE$I)←I+1)))
[NIL (if CODES←(GROUP.HIDDENCODES A)
then (ROWBLT CBROW CODES)
(ROWBLT SHAPE (perform CODEBOOK.LENGTH CODES))
else (NOCBCOUNT←NVARS)
((SHAPE$NVARS)←STARTPTR←1)
(for I to NVARS-1 do ((SHAPE$I)←I+1]
(SHOULDNT)))
[if NOCBCOUNT gt 0
then (* there is a column without a codebook)
(DPROG ((VA [if NOCBCOUNT=NVARS
then A
else (FSELECT A (ROWPTROF 'ALL
(VFROMR (for I from 1
bind PTR←STARTPTR
(NOCBSEL ←(create ROWINT
(
NELTS←NOCBCOUNT)))
repeatuntil PTR←(SHAPE$PTR)
=STARTPTR
do (NOCBSEL$I←PTR)
finally (RETURN NOCBSEL] ARRAY))
(* The non-codebook observations are read by selecting an array containing just those columns.
The codebook is constructed in CAR of CBROW. By going around the circular threaded list in SHAPE, the right position
in CBROW is found.)
[bind E (GSBVA ←(SETUP VA 'ROWMAJOR))
declare (E SCALAR)
(GSBVA GENSTATEBLOCK)
for PTR←STARTPTR by SHAPE$PTR until GSBVA:DONE
when (AND E←(GETAELT VA (NEXT GSBVA))
~(perform CODEBOOK.FINDLAB (CBROW$PTR)
E))
do (push (CBROW$PTR)
(create CODEPAIR
CODE ←(if (type? FLOATING E)
then (FPLUS E)
else (IPLUS E))
CODELAB ←(PACK* E '=])
(bind PTR←STARTPTR repeatuntil PTR=STARTPTR
do
(* Define an empty object if there were no observations for one of the none-codebook classifiers.
Sort the code book so levels will appear in numerical order.)
(PTR←(PROG1 SHAPE$PTR (SHAPE$PTR)←(perform CODEBOOK.LENGTH
(CBROW$PTR←(SORT CBROW$PTR T])])
(MOMENTS
[ULAMBDA ((A (EXPECTS ARRAY))
(WT (EXPECTS ARRAY (ONEOF NIL ARRAY)))
(M [ONEOF NIL (INTEGER (SATISFIES ~(MINUSP M]))
(* rmk: "13-JUL-78 15:47" posted: "13-JUL-78 15:47")
(* The user entry for MOMENTS. Dispatches to MOMENTS.FEW for less than 3 moments with or without weighting, to
MOMENTS.MANY for higher order moments without weighting)
(if M=NIL
then M←2)
(DPROG ((MOMEN NIL VECTOR)
(RETURNS SIMARRAY))
[MOMEN←(VFROMR (if (VSCALARP A)
then (DPROG ((MR NIL ROWSCALAR)
(N (if WT
then (CONV.SCALAR WT)
else 1.0) SCALAR)
(V (GETAELT A (VSCALARPTR A)) SCALAR)
(RETURNS ROWSCALAR))
(if (AND V N ~(MINUSP N))
else (N←0.0)
(V←NIL))
(MR←(create ROWFLOAT
NELTS ←(M+1)
INIT ←(AND V 0.0)))
(MR$1←N)
(if M GT 1
then (MR$2←V))
(RETURN MR))
else (if M GT 2
then (if WT
then (UERROR
"Weighting not implemented for higher-order moments"))
(MOMENTS.MANY A M)
else (MOMENTS.FEW A WT M]
[if LABPROPFLAG
then (SETTITLE MOMEN (MAKETITLE NIL A (if WT
then < ", Weighted by " WT>)))
(SETDIMLAB MOMEN 1 'Moment)
(for LEV to M+1 do (SETLEVLAB MOMEN 1 LEV (SELECTQ LEV
(1 'N)
(2 'Mean)
(3 'Variance)
(4 '3rd-mnt)
(PACK* LEV-1 'th-mnt]
(RETURN MOMEN))])
(MOMENTS.FEW
[DLAMBDA ((A ARRAY)
(WT (ONEOF NIL ARRAY))
(M INTEGER (SATISFIES (BETWEEN M 0 2)))
(RETURNS ROWSCALAR))
(* jop: " 5-Sep-85 17:12" posted: "25-MAY-77 16:39")
(AND WT [NOT (IEQP (NLOGICALELTS (fetch SHAPE of WT))
(NLOGICALELTS (fetch SHAPE of A]
(UERROR "Wrong size weighting array"))
(bind VAL DEV OLDNCELL (WVAL ← 1.0)
(MEAN ← 0.0)
(SS ← 0.0)
(NCELL ← 0.0)
(GSBA ←(SETUP A (QUOTE ROWMAJOR)))
[GSBW ←(AND WT (SETUP WT (QUOTE ROWMAJOR]
(MROW ←(create ROWFLOAT
NELTS ←(PLUS M 1)
INIT ← 0.0))
declare (WVAL SCALAR)
(VAL SCALAR)
(DEV FLOATING (* Deviation of current VAL from previous mean)
)
(MEAN FLOATING)
(SS FLOATING)
(GSBA GENSTATEBLOCK)
(GSBW (ONEOF NIL GENSTATEBLOCK))
(NCELL FLOATING)
(OLDNCELL FLOATING)
(MROW ROWSCALAR)
(RETURNS ROWSCALAR)
until (fetch DONE of GSBA)
do (if (NULL WT)
elseif (AND (SETQ WVAL (GETAELT WT (NEXT GSBW)))
(FGREATERP (SETQ WVAL (FLOAT WVAL))
0.0))
else (SKIP GSBA 1)
(GO $$ITERATE))
(if [NOT (SETQ VAL (GETAELT A (NEXT GSBA]
then (GO $$ITERATE))
(SETQ OLDNCELL NCELL)
(fadd NCELL WVAL)
(if (EQ M 0)
then (GO $$ITERATE))
(ASSERT (GREATERP NCELL 0.0))
(SETQ VAL (FLOAT VAL))
(DPROGN ((VAL FLOATING))
(fadd MEAN (FQUOTIENT (FTIMES WVAL (SETQ DEV (FDIFFERENCE VAL MEAN)))
NCELL))
(if (EQ M 1)
then (GO $$ITERATE))
(fadd SS (FQUOTIENT (FTIMES WVAL OLDNCELL DEV DEV)
NCELL)))
finally (SETRELT MROW 1 NCELL)
(if (EQ M 0)
elseif (GREATERP NCELL 0.0)
then (SETRELT MROW 2 MEAN)
[if (AND (EQ M 2)
(FGREATERP NCELL 1.0))
then (SETRELT MROW 3 (FQUOTIENT SS (DIFFERENCE NCELL 1.0]
else (SETRELT MROW 2 NIL)
(SETRELT MROW 3 NIL)) (* Convert sum of squares to unbiased variance
estimate)
(RETURN MROW))])
(MOMENTS.MANY
[DLAMBDA ((A ARRAY)
(M INTEGER (SATISFIES (M gt 1)) (* Number of moments-1.))
(RETURNS ROWSCALAR))
(* rmk: " 2-APR-79 17:33" posted: " 9-JUN-77 13:44")
(DPROG ((M1 (M+1) INTEGER)
THEN (OUTPUT (create ROWFLOAT
NELTS ← M1
INIT ← 0.0) ROWSCALAR)) (* binomial coefficients)
[until GSBA:DONE bind N V Z Y (BIN ←(create ROWFLOAT (NELTS←M1*(M1+1)/2)
(INIT←1.0)))
(SAVROW ←(create ROWFLOAT (NELTS←M1)
(INIT←0.0)))
(YPOW ←(create ROWFLOAT (NELTS←M1)
(INIT←1.0)))
(GSBA ←(SETUP A 'ROWMAJOR))
declare (OUTPUT ROWFLOAT (* Constrain the type))
(N FLOATING (* Number of values))
(V SCALAR (* Current value))
(Y FLOATING)
(GSBA GENSTATEBLOCK)
(Z FLOATING)
(BIN ROWFLOAT (* Binomial coefficients))
(SAVROW ROWFLOAT)
(YPOW ROWFLOAT)
first [for I from 2 to M as K from 4 by 2 do (for J from 2 to I
do (add K 1)
(BIN$K←(FPLUS BIN$(K-I)
BIN$(K-I-1]
do (AND (until V←(GETAELT A (NEXT GSBA)) do (if GSBA:DONE
then (RETURN T)))
(RETURN))
(swap SAVROW OUTPUT) (* Continue to the next non-NIL value)
(N←SAVROW$1)
((OUTPUT$1)←Y←N+1.0)
(Z←SAVROW$2)
(Y←(FQUOTIENT (FDIFFERENCE (FLOAT V)
Z)
Y)) (* y= x-m/ n)
((OUTPUT$2)←(FPLUS Y Z))
(Z←(FTIMES N Y))
(YPOW$2←Y)
(for K X B U SGN←-1.0 from 2 to M
declare (X FLOATING)
(B INTEGER)
(K INTEGER)
(SGN (MEMQ -1.0 1.0))
(U (MEMQ -1.0 1.0))
do (YPOW$(K+1)←(FTIMES YPOW$K Y))
(X←Z←(FTIMES N Y Z))
(B←(K+1)*(K+2)/2)
(U←SGN←(FMINUS SGN))
(for L from K+1 by -1 to 1 as P from 1
declare (L INTEGER)
(P INTEGER)
do (if L~=K
then (fadd X (FTIMES SAVROW$P YPOW$L BIN$B U)))
(U←(FMINUS U))
(add B -1)
finally ((OUTPUT$(P-1))←X]
(DPROG ((N (OUTPUT$1) FLOATING))
(if N gt 0.0
then (for J from 4 to M1 declare (OUTPUT ROWFLOAT)
first ((OUTPUT$3)←(if N gt 1.0
then (FQUOTIENT OUTPUT$3 N-1.0)
else 0.0))
do ((OUTPUT$J)←(FQUOTIENT OUTPUT$J N)))
(* Scaling by N for third and higher moments, N-1 for
the variance)
else (for J from 2 to M1 do ((OUTPUT$J)←NIL))))
(* If no observations, make all moments be NIL)
(RETURN OUTPUT))])
(PAIRN
[ULAMBDA ((A (EXPECTS MATRIX))
(WT (EXPECTS VECTOR (ONEOF NIL VECTOR)) (* Weight vector))
(RETURNS MATRIX))
(* edited: "16-Feb-86 12:38" posted: "28-MAY-77 00:05"
)
(if (AND WT (NEQ (GETRELT (fetch SHAPE of WT)
1)
(GETRELT (fetch SHAPE of A)
1)))
then (UERROR "Wrong length weighting vector")) (* Returns the pairwise N corresponding to a covar
matrix on A.)
(DPROG ((PAIRN NIL SIMARRAY (* The output matrix))
(NVARS (GETRELT (fetch SHAPE of A)
2) INTEGER (* Number of variables))
(N (if (NULL WT)
then 0
else (ZEROFORARRAY WT)) ARITH (* Sum of the weights))
(NROW NIL (ONEOF ROWINT ROWFLOAT) (* If needed, row of pairwise N's))
THEN (NR (IQUOTIENT (ITIMES NVARS (ADD1 NVARS))
2) IJK))
[bind MFLAG OBS (WVAL ← 1)
(GSBA ←(SETUP A (QUOTE ROWMAJOR)))
[GSBW ←(AND WT (SETUP WT (QUOTE ROWMAJOR]
declare (MFLAG BOOL)
(OBS ROWPTR (* If needed, elements=NIL if var missing)
)
(WVAL SCALAR)
(GSBA GENSTATEBLOCK)
(GSBW (ONEOF NIL GENSTATEBLOCK))
until (fetch DONE of GSBA)
do (if (NULL WT)
elseif (AND (SETQ WVAL (GETAELT WT (NEXT GSBW)))
(GREATERP WVAL 0))
else (SKIP GSBA NVARS)
(GO $$ITERATE)) (* Read observations for each var)
(for J to NVARS
do (if (NULL (GETAELT A (NEXT GSBA)))
then (if (NULL MFLAG)
then (SETQ MFLAG T)
(* The very first missing observation)
(SETQ OBS (create ROWPTR
NELTS ← NVARS
INIT ← T))
(SETQ NROW (create ROWSCALAR
NELTS ← NR
RELTTYPE ←(if (type?
FLOATING
WVAL)
then
(QUOTE
FLOATING)
else
(QUOTE INTEGER))
INIT ← N)))
(SETRELT OBS J NIL)
elseif MFLAG
then (SETRELT OBS J T)))
(if (NULL MFLAG)
then (SETQ N (PLUS N WVAL))
else (for J to NVARS bind (JK ← 0) declare (JK IJK)
do (if (GETRELT OBS J)
then [for K to J
do (add JK 1)
(if (GETRELT OBS K)
then (SETRELT NROW JK
(PLUS WVAL
(GETRELT NROW JK]
else (add JK J]
[SETQ PAIRN (create SIMARRAY
FORMAT ←(QUOTE SYMMETRIC)
SHAPE ←(ROWINTOF NVARS NVARS)
ELEMENTBLOCK ←(OR NROW
(if (type? FLOATING N)
then (create ROWFLOAT
NELTS ← NR
INIT ← N)
else (create ROWINT
NELTS ← NR
INIT ← N]
(* We're just smashing NROW into the array frame,
assuming that we computed a size consistent with the
shape)
(if LABPROPFLAG
then [SETTITLE PAIRN (MAKETITLE NIL A (if WT
then (LIST ", Weighted by " WT]
(LAB.COPYDIM A PAIRN 2 1)
(LAB.COPYDIM PAIRN PAIRN 1 2) (* Guarantee label sharing))
(RETURN PAIRN))])
(POOL
[ULAMBDA ((MTABLE (EXPECTS ARRAY) (SATISFIES MTABLE:FORMAT='FULL MTABLE:SHAPE$(MTABLE:NDIMS)=3)
"Invalid moments table")
(RETURNS VECTOR))
(* jop: " 5-Sep-85 17:19" posted: "23-SEP-77 14:11")
(* Collapses a moments table into a 3-vector of the pooled cell N, cell Means, and mean-centered sum of squares.
Effectively removes factors from a moments table.)
(bind C D TEMP NEWN POOL (LASTN ← 0.0)
(LASTMEAN ← 0.0)
(LASTSS ← 0.0)
(GSBMOM ←(SETUP MTABLE (QUOTE ROWMAJOR)))
declare (LASTN (ONEOF NIL FLOATING))
(LASTMEAN (ONEOF NIL FLOATING))
(LASTSS (ONEOF NIL FLOATING))
(C FLOATING)
(D FLOATING)
(TEMP SCALAR)
(NEWN SCALAR)
(GSBMOM GENSTATEBLOCK)
(POOL VECTOR)
until (fetch DONE of GSBMOM)
do (if [NULL (SETQ NEWN (GETAELT MTABLE (NEXT GSBMOM]
then (SETQ LASTN (SETQ LASTMEAN (SETQ LASTSS NIL)))
(RETURN)
elseif (FGREATERP NEWN 0.0)
then (DPROGN ((LASTN FLOATING)
(NEWN ARITH))
(SETQ C (FTIMES LASTN NEWN))
[if (NULL LASTMEAN)
then (SKIP GSBMOM 2)
elseif [NULL (SETQ TEMP (GETAELT MTABLE (NEXT GSBMOM]
then (SETQ LASTMEAN (SETQ LASTSS NIL))
(SKIP GSBMOM 1)
else (DPROGN ((LASTMEAN FLOATING))
(fadd LASTN NEWN)
(SETQ D (FDIFFERENCE TEMP LASTMEAN))
(if [AND LASTSS (SETQ TEMP (GETAELT MTABLE (NEXT GSBMOM]
then (DPROGN ((LASTSS FLOATING)
(TEMP ARITH))
[fadd LASTSS (FTIMES TEMP (DIFFERENCE NEWN 1.0))
(FTIMES C D (SETQ D (FQUOTIENT D LASTN])
else (SETQ LASTSS NIL)
(SETQ D (FQUOTIENT D LASTN)))
(fadd LASTMEAN (FTIMES NEWN D)))])
else (SKIP GSBMOM 2))
finally [SETQ POOL (VFROMR (FILLROW (create ROWSCALAR
NELTS ← 3
RELTTYPE ←(QUOTE FLOATING))
1 LASTN LASTMEAN (AND LASTN LASTSS
(FQUOTIENT LASTSS (DIFFERENCE LASTN
1.0]
(* If LASTN=1, then LASTSS will be zero, and the LISP
FQUOTIENT will return 0.0)
(if LABPROPFLAG
then (LAB.COPYDIM MTABLE POOL (fetch NDIMS of MTABLE)
1)
(SETTITLE POOL (MAKETITLE NIL MTABLE)))
(RETURN POOL))])
)
(PUTPROPS COMPRESS COPYRIGHT ("Xerox Corporation" 1983 1984 1985 1986))
(DECLARE: DONTCOPY
(FILEMAP (NIL (575 41758 (COUNTS 585 . 1413) (COVAR 1415 . 11819) (GROUP 11821 . 19798) (
GROUP.HIDDENCODES 19800 . 20499) (GROUP.HIDDENLEVLAB 20501 . 21413) (GROUP.LAB 21415 . 22492) (
GROUP.LINLOC 22494 . 23456) (GROUP.SHAPE 23458 . 27401) (MOMENTS 27403 . 29400) (MOMENTS.FEW 29402 .
31865) (MOMENTS.MANY 31867 . 35084) (PAIRN 35086 . 38926) (POOL 38928 . 41756)))))
STOP