Change row values to zero if less than row standard deviation
Asked Answered
P

2

5

I want to change all values of a row to zero if they are less than the standard deviation of that row.

set.seed(007)
X <- data.frame(matrix(sample(c(5:50), 100, replace=TRUE), ncol=10))

   X1 X2 X3 X4 X5 X6 X7 X8 X9 X10
1  37 10 43 45 11 17 39 13 13  44
2  10 24 32 16  7 50 41 47  9  39
3  23 49 46 35 16 30 22 10 11  46
4  41 46 19 28 47 39 27 40 49  13
5  29 23 49 10 50 17 42 43  7  31
6  31 26 11 36 35 43 45 29 33   9
7  21 12  5 21 29 12 31 30  7  30
8  32 24  8 43  9 17 35 44 41   8
9  20 44 39  8 40 17 27 45 14  37
10 50  8  5 48 27 15 15 12 30  15

The lines below appear to do the job, but is terribly slow on my actual use-case and I'm a bit unsure what sapply is returning....

Y <- t(sapply(1:nrow(X), function(i) 
      sapply(1:ncol(X), function(j) 
        ifelse(X[i,][[j]] < sd(X[i,]), 0, X[i,][[j]]))))

      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]   37    0   43   45    0   17   39    0    0    44
 [2,]    0   24   32    0    0   50   41   47    0    39
 [3,]   23   49   46   35   16   30   22    0    0    46
 [4,]   41   46   19   28   47   39   27   40   49    13
 [5,]   29   23   49    0   50   17   42   43    0    31
 [6,]   31   26    0   36   35   43   45   29   33     0
 [7,]   21   12    0   21   29   12   31   30    0    30
 [8,]   32   24    0   43    0   17   35   44   41     0
 [9,]   20   44   39    0   40   17   27   45   14    37
[10,]   50    0    0   48   27    0    0    0   30     0

What is a faster and more efficient method?

UPDATE Thank you all very much for the quick and efficient answers!

Here's how they stack up...

set.seed(007)
size <- 1e5
X <- matrix(sample(c(5:50), size, replace=TRUE), ncol=size/2)

library(microbenchmark)    
results <- microbenchmark(
  X[ sweep(X, 1, apply(X,1,sd) ) < 0 ] <- 0,
  X[t(apply(X, 1, function(x) x - sd(x) < 0))] <- 0,
  sapply(X, function(x) ifelse(x < sd(x), 0, x)),
  times = 100L)
print(results)
Unit: milliseconds
                                              expr         min          lq     median         uq        max neval
          X[sweep(X, 1, apply(X, 1, sd)) < 0] <- 0    7.966167   10.869785   12.38399   15.00107   45.41557   100
 X[t(apply(X, 1, function(x) x - sd(x) < 0))] <- 0    7.344227    9.675577   11.22283   14.34280   53.70728   100
    sapply(X, function(x) ifelse(x < sd(x), 0, x)) 3028.336236 3221.325598 3302.16115 3466.66875 4539.88358   100
# plot
if (require("ggplot2")) {
  plt <- ggplot2::qplot(y=time, data=results, colour=expr)
  plt <- plt + ggplot2::scale_y_log10()
  print(plt)
}

enter image description here

Looks like Arun's answer is the fastest by a tiny bit (as Arun notes). However, DWin's is eight characters less typing and is notable for using the exotic (to me) sweep function.

A minor recreational update, Arun's method is significantly faster (t = 2.0112, df = 191.985, p-value = 0.04571) or, if you prefer, the mean speed of Arun's function is credibly faster than the mean speed of DWin's (using this robust Bayesian estimation method, Group 1 = DWin, Group 2 = Arun, though Arun's timings are not a good fit for the t-dist):

enter image description here

Provision answered 16/4, 2013 at 20:50 Comment(3)
(+1) for the seed number.. :)Sabbath
:) taken from https://mcmap.net/q/267064/-calculating-standard-deviation-of-each-rowProvision
At one time the discovery of the sweep operation was a big deal, but I think its novelty has faded.Hungary
H
3

I suspect this is slower that the apply solution, but since there is no need to add the data.frame step and the fact that apply.data.frame is notoriously slow, I may still "win" or "keep even" at least until the other contestants tumble to the fact that I use a matrix object.

set.seed(007)
X <- matrix(sample(c(5:50), 100, replace=TRUE), ncol=10)
X[ sweep(X, 1, apply(X,1,sd) ) < 0 ] <- 0

Note that Richardo and I both got the same different starting point than the OP although I think he needed to transpose if he wants a row operation:

> X
   X1 X2 X3 X4 X5 X6 X7 X8 X9 X10
1  50  0 34 36 41 31  0 18 45  20
2  23 15 18 17 22 38 28 32 45   0
3   0 40 50  0 39 40 40 43 16  46
4   0  0 46  0 25 33 36 33 39   0
5  16 25 50 22 46 38 30  0 22  38
6  41  0  0 43 19 22 35 31  0  31
7  20 30 33 27  0 12 26 25  0  29
8  49  0 27 41 42  0 27 25 40  21
9   0 50 49 43 46 22 20 33 21  42
10 26 19 21 26 49 17 24 47 24  13

Added note: I was playing around with the rowMeans function to see if I could come up with a vectorized alternative to apply(X,1,sd) version of sd():

sqrt(rowSums((X[1:10, ]-rowMeans(X))^2)/9)

So:

 sdbyrow <- function(mat) sqrt(rowSums((mat-rowMeans(mat))^2)/(ncol(mat)-1) )
 all.equal(apply(X,1,sd), sdbyrow(X) )
#[1] TRUE
Hungary answered 16/4, 2013 at 21:12 Comment(4)
The apply version with the transpose (t) seems to edge sweep a bit in my benchmarking on a 1e5 by 1e3 data.frame. (2.37 vs 2.72 secs)Sabbath
If I may be permitted a minor further clarification, if I wanted to keep only the row values that are less than one standard deviation above the minimum value of the row, would this be correct: XX[ sweep(XX, 1, (apply(XX,1,min) + apply(XX,1,sd)) ) < 0 ] <- 0?Provision
It looks correct to me. There is also a rowMin in one of the contributed packages that might do the same as apply(X,1,min) and perhaps even a rowSD, although I haven't checked that possibility yet.Hungary
Great, thanks. I've upvoted a bunch of your other recent answers as a token of gratitude since upvoting comments is worth nothing (done it anyway)Provision
S
4

How about this?

X[t(apply(X, 1, function(x) x - sd(x) < 0))] <- 0
#    X1 X2 X3 X4 X5 X6 X7 X8 X9 X10
# 1  50  0 34 36 41 31  0 18 45  20
# 2  23 15 18 17 22 38 28 32 45   0
# 3   0 40 50  0 39 40 40 43 16  46
# 4   0  0 46  0 25 33 36 33 39   0
# 5  16 25 50 22 46 38 30  0 22  38
# 6  41  0  0 43 19 22 35 31  0  31
# 7  20 30 33 27  0 12 26 25  0  29
# 8  49  0 27 41 42  0 27 25 40  21
# 9   0 50 49 43 46 22 20 33 21  42
# 10 26 19 21 26 49 17 24 47 24  13
Sabbath answered 16/4, 2013 at 20:54 Comment(0)
H
3

I suspect this is slower that the apply solution, but since there is no need to add the data.frame step and the fact that apply.data.frame is notoriously slow, I may still "win" or "keep even" at least until the other contestants tumble to the fact that I use a matrix object.

set.seed(007)
X <- matrix(sample(c(5:50), 100, replace=TRUE), ncol=10)
X[ sweep(X, 1, apply(X,1,sd) ) < 0 ] <- 0

Note that Richardo and I both got the same different starting point than the OP although I think he needed to transpose if he wants a row operation:

> X
   X1 X2 X3 X4 X5 X6 X7 X8 X9 X10
1  50  0 34 36 41 31  0 18 45  20
2  23 15 18 17 22 38 28 32 45   0
3   0 40 50  0 39 40 40 43 16  46
4   0  0 46  0 25 33 36 33 39   0
5  16 25 50 22 46 38 30  0 22  38
6  41  0  0 43 19 22 35 31  0  31
7  20 30 33 27  0 12 26 25  0  29
8  49  0 27 41 42  0 27 25 40  21
9   0 50 49 43 46 22 20 33 21  42
10 26 19 21 26 49 17 24 47 24  13

Added note: I was playing around with the rowMeans function to see if I could come up with a vectorized alternative to apply(X,1,sd) version of sd():

sqrt(rowSums((X[1:10, ]-rowMeans(X))^2)/9)

So:

 sdbyrow <- function(mat) sqrt(rowSums((mat-rowMeans(mat))^2)/(ncol(mat)-1) )
 all.equal(apply(X,1,sd), sdbyrow(X) )
#[1] TRUE
Hungary answered 16/4, 2013 at 21:12 Comment(4)
The apply version with the transpose (t) seems to edge sweep a bit in my benchmarking on a 1e5 by 1e3 data.frame. (2.37 vs 2.72 secs)Sabbath
If I may be permitted a minor further clarification, if I wanted to keep only the row values that are less than one standard deviation above the minimum value of the row, would this be correct: XX[ sweep(XX, 1, (apply(XX,1,min) + apply(XX,1,sd)) ) < 0 ] <- 0?Provision
It looks correct to me. There is also a rowMin in one of the contributed packages that might do the same as apply(X,1,min) and perhaps even a rowSD, although I haven't checked that possibility yet.Hungary
Great, thanks. I've upvoted a bunch of your other recent answers as a token of gratitude since upvoting comments is worth nothing (done it anyway)Provision

© 2022 - 2024 — McMap. All rights reserved.