Cross Tabulations
Vector 19:1, page 29
Transcribed from the printed issue; not yet reviewed. Read it in the PDF of the issue, page 29
Utility Corner: Cross Tabulations¶
(The purpose of this column is to make you more productive by introducing you to utility functions. Think of utility functions as APL functions that have names instead of symbols. By expanding your function vocabulary, you’ll be able to write APL code that’s more concise, more efficient, and more readable.)
Suppose you have a "database" of employee information. For each of the 2000 employees, you have four items of information: age, state of residence, type of employee, and bonus amount received last year. Here are a few records:
The information is stored in your APL workspace in four simple variables: AGE, STATE, TYPE, and BONUS. AGE and BONUS are 2000-element numeric vectors; STATE is a 2000-row, 2-column character matrix; and TYPE is a 2000-element character vector.
In last issue’s Limbering Up column, we asked you to summarise this mass of information by grouping the employees by age, state, and type. We gave you the following function to complete:
∇ R←CROSSTAB
[1] ⍝ Summarises employee information
[2] ⍝ stored in workspace globals AGE,
[3] ⍝ STATE, TYPE, and BONUS. AGE and
[4] ⍝ BONUS are numeric vectors, STATE
[5] ⍝ is a two-column character matrix,
[6] ⍝ and TYPE is a character vector.
[7] ⍝ Their elements/rows are in 1-to1
[8] ⍝ correspondence. R is a four-
[9] ⍝ dimensional array whose shape is
[10] ⍝ (5 10 3 2):
[11] ⍝
[12] ⍝ R[I;;;] The Ith age range:
[13] ⍝ I=1; 0≤AGE≤29
[14] ⍝ I=2: 30≤AGE≤39
[15] ⍝ I=3: 40≤AGE≤49
[16] ⍝ I=4: 50≤AGE≤64
[17] ⍝ I=5: 65≤AGE
[18] ⍝ R[;J;;] The Jth state. The glo-
[19] ⍝ bal ALLSTATES is a 10-row,
[20] ⍝ 2-column character matrix of
[21] ⍝ the codes for all ten states
[22] ⍝ in which the employees reside.
[23] ⍝ R[;;K;] The type 'EOS'[K]
[24] ⍝ R[;;;1] The employee count.
[25] ⍝ R[;;;2] The total bonus amount.
∇
The idea behind a function like CROSSTAB is to process the database once, summarising its information into a single multidimensional array. All subsequent analysis can then be efficiently performed on the array, rather than reprocessing the thousands (or millions) of records.
The typical APL approach to solving multidimensional cross-tabulation problems like the one above is to use outer and inner products. For example, you can use ∘.= to get the number of employees by type:
It’s a small step to get the bonus amount by type:
By laminating a row of 1s onto BONUS, you can get both employee count and bonus amount with one expression:
Likewise, you can get employee count and bonus amount by state:
How do we put the two together to get employee count and bonus amount by state? Unfortunately, the inner/outer product approach gets much harder to formulate as you go beyond two dimensions. Here’s a painful, but successful, attempt:
This expression requires the reader to do a lot of head-scratching; and it requires the computer to do a lot of work. Instead of getting out the crow-bar to pry the final dimension (age) into the expression, let’s look for a different approach.
One possibility is to use nesting:
(Replace the APL*PLUS split ↓ and mix ↑ function by ⊂[2] and ⊃[2] if using APL2.)
It is debatable whether this expression requires more or less head-scratching. Moreover, it is not more efficient than the previous expressions. All the same comparisons are being made. However, it does lend itself to be extended to higher dimensions. (Bet you can hardly wait!)
This expression solves our problem, but is extremely inefficient. Let’s formulate a different approach by working backwards from the desired result. Where do each of the original records affect the eventual result?
AGE/ STATE/ TYPE/ Result
ind. ind. ind. ind. BONUS
26/1 FL/5 E/1 13 300
39/2 MI/8 S/3 54 550
31/2 NY/3 O/2 38 720
52/4 NY/3 S/3 99 375
In this table, we’ve redisplayed the four sample records from above. New each AGE, STATE, and TYPE value, we’ve included the index of the class implied by that value. For example, look at the last record. Age 52 is in the 4th age class (50 to 64); NY is in the 3rd state class (i.e. ALLSTATES[3;]='NY'); and S is the 3rd employee type (E=1, O=2, S=3).
In an array of 5 (age) by 10 (state) by 3 (type) possibilities, the cell affected this last record is the one in the 4th plane, 3rd row, and 3rd column, or the 99th cell out of 150:
By translating the age, state, and type values into class indices, and then combining them into a single "result index", the multi-dimensional problem has been reduced into a one-dimensional problem. You have all the information you need to solve the problem in the last two columns of the table above.
Accumulate the BONUS amounts and record counts by each distinct result index; assign the accumulated amounts and counts into their corresponding positions of two 150-element vectors (or a 150 by 2 matrix); and reshape the result into the desired 4 dimensions.
The accumulation logic is the standard (and efficient) sort, shift, and compare logic that’s been around for eons. Here’s the finished code:
∇ R←CROSSTAB;AI;B;C;G;L;RI;SI;SRI;TI
[1] ⍝ Summarises employee information
: : : :
[26] ⍝ AGE class index:
[27] AI←0 30 40 50 65 LIOTAI AGE
[28] ⍝ STATE class index:
[29] SI←ALLSTATES CMIOTA STATE
[30] ⍝ TYPE class index:
[31] TI←'EOS'⍳TYPE
[32] ⍝ Result index:
[33] RI←(30×AI-1)+(3×SI-1)+TI
[34] ⍝ Sorted result indices:
[35] SRI←RI[G←⍋RI]
[36] ⍝ Flag last of each like index:
[37] L←SRI≠1↓SRI,¯1
[38] ⍝ Distinct result indices:
[39] SRI←L/SRI
[40] ⍝ Cumulative employee counts
[41] C←L/⍳⍴L
[42] ⍝ Employee counts
[43] C←C-¯1↓0,C
[44] ⍝ Cumulative bonus amounts:
[45] B←L/+\BONUS[G]
[46] ⍝ Bonus amounts:
[47] B←B-¯1↓0,B
[48] ⍝ Construct result; insert; reshape:
[49] R←150 2⍴0
[50] R[SRI;]←C,[1.5]B
[51] R←5 10 3 2⍴R
∇
Notice the use of LIOTAI (lower limit iota, inclusive), presented in the previous newsletter, and the of CMIOTA (character matrix iota), presented several issues back, to translate the age and state values to indices.
The logic in the above CROSSTAB function has been generalised and included in a utility function named PLUSRED. Here’s how you’d write CROSSTAB using PLUSRED:
∇ R←CROSSTAB;AI;SI;TI
[1] ⍝ Summarises employee information
: : : :
[26] ⍝ AGE class index:
[27] AI←0 30 40 50 65 LIOTAI AGE
[28] ⍝ STATE class index:
[29] SI←ALLSTATES CMIOTA STATE
[30] ⍝ TYPE class index:
[31] TI←'EOS'⍳TYPE
[32] ⍝ The data to summarise (incl freq):
[33] R←1,[0.5]BONUS
[34] R←5 10 3 AI SI TI PLUSRED R
∇
The left argument of PLUSRED is a nested argument. (A non-nested version of PLUSRED can be found in the Zark Library of Utility Functions.) The first three items are the class sizes. There are 5 age classes, 10 state classes, and 3 type classes. The next three items are the respective vectors of class indices. The right argument is the array to be summarised; its last dimension is "reduced," i.e. replaced by the class sizes. In this case, if the argument has shape 2 2000, the result has shape 2 5 10 3. If the shape of the argument was instead 2000 2, we would want to "reduce" the first dimension (the 2000), and get back a result with shape 5 10 3 2. To specify a reduction dimension other than the last, include its index at the end of PLUSRED’s left argument, like so:
Here’s the PLUSRED utility function:
∇ R←L PLUSRED ARRAY;CUM;DIM;DSHAPE;G;I;LAST;M;N;RANK;RRI;S;U
[1] ⍝ No. of ways for N-way reduction:
[2] N←⌊(×/⍴L)÷2
[3] ⍝ Which dimension to be "reduced"?
[4] DIM←1↑((N+N)↓L),¯1↑⍳RANK←⍴⍴ARRAY
[5] ⍝ Get replacement shape from left arg:
[6] DSHAPE←N⍴L
[7] ⍝ Initial "raveled reslt inds":
[8] I←⎕IO
[9] ⍝ Index from ⍳DSHAPE[I] to cause
[10] ⍝ index err if invalid indices:
[11] RRI←(⍳DSHAPE[I])[(N+I)⊃L]
[12] ⍝ Continue computing RRI (N loops
[13] ⍝ for N-way reduction):
[14] L1:→((N+⎕IO)≤I←I+1)⍴L2
[15] RRI←(⍳DSHAPE[I])[(N+I)⊃L]+DSHAPE[I]×RRI-⎕IO
[16] →L1
[17] ⍝ Determine unique elements of RRI:
[18] L2:S←RRI[G←⍋RRI]
[19] U←(LAST←S≠1↓S,¯1)/S
[20] ⍝ Reorder ARRAY to conform with S:
[21] M←DIM-⎕IO
[22] ARRAY←⍎'ARRAY[',(M⍴';'),'G',((RANK-M+1)⍴';'),']'
[23] ⍝ Perform partitioned reduction:
[24] CUM←LAST/[DIM]+\[DIM]ARRAY
[25] CUM←CUM-(⍴CUM)↑0,[DIM]CUM
[26] ⍝ Initialise result. Fill with 0s
[27] ⍝ Ravel the DSHAPE dim.s:
[28] R←((-M)⌽(×/DSHAPE),1↓M⌽⍴ARRAY)⍴0
[29] ⍝ Insert result of part. reduction:
[30] ⍎'R[',(M⍴';'),'U',((RANK-M+1)⍴';'),']←CUM'
[31] ⍝ Reshape to desired shape:
[32] R←((-M)⌽DSHAPE,1↓M⌽⍴ARRAY)⍴R
∇
In case you missed the last few issues, here are the LIOTAI and CMIOTA utility functions:
∇ R←LOWER LIOTAI VALS;G;I;M
[1] ⍝ Returns the index of the value of
[2] ⍝ the vector LOWER that is the
[3] ⍝ nearest lower bound (inclusive)
[4] ⍝ for each element of the vector
[5] ⍝ VALS. Values in LOWER need not be
[6] ⍝ ascending or distinct. ⎕IO+⍴LOWER
[7] ⍝ is returned for elements in VALS
[8] ⍝ above the smallest value in LOWER
[9] ⍝
[10] ⍝ Combine and sort up to use 1st val
[11] ⍝ in LOWER if vals not distinct:
[12] R←G←⍒VALS,LOWER
[13] ⍝ Flag elts from LOWER in grade vec:
[14] I←(G≥⎕IO+⍴VALS)/⍳⍴G
[15] ⍝ Inds into LOWER, in sorted order,
[16] ⍝ using "one greater" for values
[17] ⍝ above largest in LOWER:
[18] M←(G[I]-⍴VALS),⎕IO+⍴LOWER
[19] ⍝ Replicate inds for values in
[20] ⍝ corresp. ranges:
[21] R[G]←((I,⎕IO+⍴G)-⎕IO,I)/M
[22] R←(-⍴LOWER)↓R ⍝ Discard LOWER inds
∇
∇ R←A CMIOTA1 B;G;J;K;N;S;T;Z
[1] ⍝ Returns row indices of character
[2] ⍝ matrix A where each row of B is
[3] ⍝ found. Like dyadic ⍳, CMIOTA re-
[4] ⍝ turns one greater than the largest
[5] ⍝ index in A for rows of B not found.
[6] ⍝
[7] ⍝ Index returnd if not found:
[8] Z←⎕IO+1↑⍴A
[9] N←1↑⍴B ⍝ No. rows in B
[10] ⍝ Stick ¯1↑⎕AV as new row of A:
[11] A←A,[⎕IO]⎕AV[255+⎕IO]
[12] G←⎕AV⍋B,[⎕IO]A ⍝ Grade vec
[13] ⍝ Which ones come from B (∧/K∊⍳N):
[14] K←(S←G<N+⎕IO)/G
[15] ⍝ The relative ind of the next A row:
[16] J←⎕IO+S/+\T←~S
[17] R←J←(T/G)[J]-N ⍝ The absolute A inds
[18] R[K]←J ⍝ ...in the orig order of B
[19] ⍝ Use the "not-found" value for those
[20] ⍝ not matching:
[21] R[(∨/B≠A[R;])/⍳N]←Z
∇