餡子付゛録゛

ソフトウェア開発ツールの便利な使い方を紹介。

Rで各集団における大きい方から2番目の値を調べる

ツイッターで見かけた御題なのですが、(学校の課題などで)ありそうなので、置いておきます。

# お試しデータを作る
set.seed(20150208)
n <- 30
df <- data.frame(type=c("A", "B", "C")[round(runif(n, min=0.5, max=3.5))], value=60 - ((1:n)-10)^2)

# タイプごとに最大値を求める
tapply(df$value, df$type, max)

# 2番目の値を求める
tapply(df$value, df$type, function(x){
  # 降順ソートをして2番目を返す
  sort(x, decreasing=TRUE)[2]
})

奥村晴彦氏から計算量がO(n log n)でオセーヨってツッコミが来たので、O(n)になるように改良してみましたが、300万件でもuserで2割ぐらいしか時間を削れませんでした*1小選挙区の数は300しかないので(謎)、遅いコードで許してください。
以下は時間計測に使ったコードです。

set.seed(20150208)

# 300万件でテスト
n <- 3000000

# データは散らばるようにしておく
df <- data.frame(type=c("A", "B", "C")[round(runif(n, min=0.5, max=3.5))], value=runif(n, min=0, max=n))

# ソートする遅い版: O(n log n)
gc();gc();system.time({
  tapply(df$value, df$type, function(x){
   sort(x, decreasing=TRUE)[2]
  })
})

# 最大値を除外する: O(n)が2回
gc();gc();system.time({
  tapply(df$value, df$type, function(x){
    max( x[x!=max(x)] )
  })
})

# コメントで示唆されたもの: O(n)が2回
gc();gc();system.time({
  tapply(df$value, df$type, function(x){
    x[which(x==max(x))] <- NA;
    max(x, na.rm=TRUE)
  })
})

*1:3000万件でも3割ぐらいしか削れないので、Rのsort関数が速く、max関数が遅いんじゃないかと言う疑惑が・・・