4

Suppose I have a data.table as below (where you can think of w as a grouping variable):

set.seed(1)
prQ = CJ(Q1 = 1:10, Q2=1:10,w=1:2)
prQ[,pQ:=runif(100,0,1)]
prQ[,pQ:=pQ/sum(pQ),by=w]

  > prQ
     Q1 Q2 w          pQ
  1:  1  1 1 0.004889560
  2:  1  1 2 0.007553012
  3:  1  2 1 0.010549565
  4:  1  2 2 0.018433927
  5:  1  3 1 0.003714138
 ---                    
196: 10  8 2 0.016183006
197: 10  9 1 0.008384253
198: 10  9 2 0.008323492
199: 10 10 1 0.014932841
200: 10 10 2 0.012278353

How do I calculate a new column, for a given w, called CDF, that does the following:

For example suppose Q1 = 4 and Q2 = 6. Define a new column,

CDF = sum(pQ) for all Q1<=4 and Q2<=6, holding w fixed.

For example, a single row:

CDF0 = sum(prQ[Q1<=4 & Q2<=6 & w==1,pQ])
prQ[Q1==4 & Q2==6,CDF:=CDF0]

I want to do this with all rows for a given w.

Desired output done using brute force:

for(w0 in 1:2){
  for(j in 1:10){
    for(p in 1:10){
          CDF0 = sum(prQ[Q1<=j & Q2<=p & w==w0,pQ])
          prQ[Q1==j & Q2==p & w==w0,CDF:=CDF0]

    }
  }
}



  > head(prQ)
   Q1 Q2 w          pQ         CDF
1:  1  1 1 0.004889560 0.004889560
2:  1  1 2 0.007553012 0.007553012
3:  1  2 1 0.010549565 0.015439125
4:  1  2 2 0.018433927 0.025986939
5:  1  3 1 0.003714138 0.019153263
6:  1  3 2 0.018234648 0.044221587
Khaynes
  • 1,976
  • 2
  • 15
  • 27
wolfsatthedoor
  • 7,163
  • 18
  • 46
  • 90
  • 1
    @PoGibas - I *think* it's more complicated than that - I think it will need to be a self join so all rows with values less the current row are summed. – thelatemail Jan 09 '19 at 23:27
  • @thelatemail Correct, PoGibas solution does not work – wolfsatthedoor Jan 09 '19 at 23:27
  • @PoGibas Ok, I have added it now. I calculated it using the brute force approach. I need something more efficient than this obviously. :) – wolfsatthedoor Jan 09 '19 at 23:39
  • I can't get it to spit out a proper result when saving to a new variable, but `prQ[ prQ, on=c("w==w","Q1<=Q1","Q2<=Q2"), sum(pQ), by=.EACHI, allow.cartesian=TRUE ]` seems close. – thelatemail Jan 10 '19 at 00:13
  • @PoGibas . It does not work for the latter rows now. If you try to compare answers using `identical(prQ1, prQ2)` where `prQ1` is the one crated with brute force and `prQ2` is the one using your code, it says `FALSE`. – wolfsatthedoor Jan 10 '19 at 00:19

1 Answers1

1

A possible approach to sum every possible sub-matrix within the matrix (with number of rows = number of unique Q2 and number of columns = number of unique Q1) constructed from pQ values:

#ensure that order is correct as values will be used to generate the matrix 
#so that all elements in the top left sub-matrix will always be 
#smaller than or equal to the bottom right element of this sub-matrix
setorder(prQ, w, Q1, Q2)

#create all possible permutations of row and column indices
subMatIdx <- prQ[, CJ(as.integer(as.factor(Q1)), as.integer(as.factor(Q2)), unique=TRUE)]

#sum every sub matrix
prQ[, CDF :=
    {
        nr <- uniqueN(Q2)

        .(Map(function(i, j) sum(matrix(pQ, nrow=nr)[1L:j, 1L:i]), 
            subMatIdx[["V1"]], subMatIdx[["V2"]]))
    },
    by=.(w)]

output:

     Q1 Q2 w          pQ        CDF
  1:  1  1 1 0.004889560 0.00488956
  2:  1  2 1 0.010549565 0.01543912
  3:  1  3 1 0.003714138 0.01915326
  4:  1  4 1 0.017396970 0.03655023
  5:  1  5 1 0.011585652 0.04813589
 ---                               
196: 10  6 2 0.001196193  0.5713282
197: 10  7 2 0.017785668  0.6535378
198: 10  8 2 0.016183006  0.7734989
199: 10  9 2 0.008323492   0.871678
200: 10 10 2 0.012278353          1

edit: what if Q1 and Q2 are negative or any real number? the line on subMatIdx should already have taken care of it.

e.g.:

set.seed(1)
prQ = CJ(Q1 = -1:10, Q2=-1:10,w=1:2)
prQ[,pQ:=runif(nrow(prQ),0,1)]
prQ[,pQ:=pQ/sum(pQ),by=w]

setorder(prQ, w, Q1, Q2)

#create all possible permutations of row and column indices
subMatIdx <- prQ[, CJ(as.integer(as.factor(Q1)), 
    as.integer(as.factor(Q2)), unique=TRUE)]

prQ[, CDF := {
        nr <- uniqueN(Q2)

        .(Map(function(i, j) sum(matrix(pQ, nrow=nr)[1L:j, 1L:i]), 
            subMatIdx[["V1"]], subMatIdx[["V2"]]))
    },
    by=.(w)]

output:

     Q1 Q2 w          pQ         CDF
  1: -1 -1 1 0.003607862 0.003607862
  2: -1  0 1 0.007784212  0.01139207
  3: -1  1 1 0.002740553  0.01413263
  4: -1  2 1 0.012836710  0.02696934
  5: -1  3 1 0.008548709  0.03551805
 ---                                
284: 10  6 2 0.011164332   0.6425251
285: 10  7 2 0.007638237   0.7360602
286: 10  8 2 0.005403923   0.8270053
287: 10  9 2 0.002008067   0.9193811
288: 10 10 2 0.002242777           1
chinsoon12
  • 25,005
  • 4
  • 25
  • 35
  • If you wanted to use a more general set of variables than w, how do you deal with this in "setorder" function? – wolfsatthedoor Jan 10 '19 at 02:22
  • you mean group by more than 1 variable? if yes, you can just add them in `by=.(w,x,y,z,...)`. the `w` in `setorder` was meant to view the original dataset for coding my next step – chinsoon12 Jan 10 '19 at 02:24
  • Like `setorder(prQ, .(w,x,y,z), Q1, Q2)`? – wolfsatthedoor Jan 10 '19 at 02:26
  • `setorder(prQ, Q1, Q2)` will do. so that all elements in the top left sub-matrix will always be smaller than or equal to the bottom right element of this sub-matrix – chinsoon12 Jan 10 '19 at 02:26
  • don't need wrap `cols` in a list (i.e. `.()`) if you have it as a char vector. just use `by=cols` – chinsoon12 Jan 10 '19 at 02:39
  • Figured it out right as you posted (thank you)...last Q, what do I need to do if Q1 and Q2 range from -1 to 10, it throws back negative subindex errors. – wolfsatthedoor Jan 10 '19 at 02:40
  • 1
    the line on `subMatIdx` should have already taken care of it. as I convert the values into integral indices. – chinsoon12 Jan 10 '19 at 02:42