ダメ出し：ちゃんと関数があるんだから，モ～

久しぶりの大物

「多次元度数分布を取る」より

プログラム自体は，やっと <- の前後に空白をおくなどで少しは読みやすくなったけど，まだ，){ や if( や for( があったり，カンマの後に空白を置かないなど相変わらず，読みにくい。

3次元クロス集計をする自前の関数を書いたというところだけど，リストや for 文を惜しげもなく使っているので，処理時間がべらぼうにかかる。

> set.seed(123456789)
> X <- matrix(rnorm(100000*3),ncol=3) + matrix(rnorm(100000*3,mean=3,sd=2),ncol=3)
>
> ns <- c(5,4,3)
> system.time(mult.hist.out <- multi.hist(X,ns = ns))
   ユーザ   システム       経過  
    21.665      0.346     21.797

これでは困る。R にはちゃんとした関数 xtabs があるのだから，それ使おうヨ。

ただ，xtabs は，元のデータがカテゴリー化されているというのが前提なので，カテゴリー化は cut 関数を使って自分で行う。

もとのプログラムも cut 関数で使う breaks と同じ働きをする変数の設定はしているが，ちょっと小細工をしている。しかし，「*0.0001」はいいかげんなものだ。意味不明のオマジナイノジュモンは使わない。

    if(is.null(bins)){
        range.X <- apply(X,2,range)
        range.X[2,] <- range.X[2,] + min(range.X[2,]-range.X[1,])*0.0001
        bins <- list()
        for(i in 1:length(X[1,])){
            bins[[i]] <- seq(from=range.X[1,i],to=range.X[2,i],length=ns[i]+1)
        }
    }

正統的なプログラムにしよう。
div 関数が元のデータのカテゴリー化を行う。cut 関数で使う breaks は，一番最初の区間の左端は含まれないので（マニュアル嫁） seq(range.x[1], range.x[2], length=ns+1) だけだと最小値の行き場所がなくなる（NA になる）ので，.Machine$double.eps を引いてやる。

mx <- function(x) {
    x <- data.frame(x)
    div <- function(x, ns) {
        range.x <- range(x)
        breaks <- seq(range.x[1], range.x[2], length=ns+1)
        breaks[1] <- breaks[1]-.Machine$double.eps
        return(as.integer(cut(x, breaks=breaks)))
    }    
    return(xtabs(~., mapply(div, x, ns)))
}

このプログラムの実行時間は

> system.time(mx(X))
   ユーザ   システム       経過  
     0.810      0.034      0.842

となり，26倍ほど速い。

自分でプログラムするにしても，もう少し速いプログラムを書ける（ようになってネ）。

ちなみに，元記事のプログラム名は multi.hist となっているが，hist はヒストグラム（度数分布図）。
度数分布表は frequency table，クロス集計表は contingency table とかいうから，名前を変えたほうがよい。