キー入力が短ければよいというものでもない

「Rで,長さ100の数列があって,連続する5個の数の和を全て(96個)欲しいとき,どうすれば賢いかな」だけど..

何をもって「賢い」とすべきかだが...

embed のソース見てみると分かるけど，結構真っ正直にやっているだけ。
だったら，ださい for を使って，もっともっと素朴にやってもおんなじようなものだろう思ったが，意外にも速かった。ベクトル演算が速いからかな。

> x = 1:10^6
> m = 5
> system.time(rowSums(embed(x, m)))
   ユーザ   システム       経過  
     0.233      0.033      0.264
> system.time({
+     n = length(x)
+     ans = x[1:(n-m+1)]
+     for (i in 2:m) {
+         ans = ans+x[i:(n-m+i)]
+     }
+ })
   ユーザ   システム       経過  
     0.074      0.012      0.085

> x = 1:10^8
> m = 5
> system.time(rowSums(embed(x, m)))
   ユーザ   システム       経過  
    15.366      3.920     19.566
> system.time({
+     n = length(x)
+     ans = x[1:(n-m+1)]
+     for (i in 2:m) {
+         ans = ans+x[i:(n-m+i)]
+     }
+ })
   ユーザ   システム       経過  
     7.871      2.709     10.390

正規表現

Rでグラフの文字をたて書きにする においてだが...

for は使っても良いが，別解を提案

tate.txt = function(x) {
  return(sapply(x, function(y) gsub("(.)", "\\1\n", y)))
}

関数の返すオブジェクトを利用しよう

メモ：独立変数にダミー変数が含まれるロジスティック回帰モデルから平均生起率を取り出すRスクリプト であるが...

「ダミー変数の標準化偏回帰係数」 で述べたように，glm 関数の返すオブジェクトを積極的に利用することにより，間違いの少ない，関数にもしやすいプログラムが書ける。

P.predi = function(object, target) {
    d = model.matrix(object$terms, eval(object$model, parent.frame()))
    x = min(d[, target]):max(d[, target])
    n = length(x)
    d2 = matrix(rep(colMeans(d), each = n), nrow = n)
    colnames(d2) = names(object$coefficients)
    d2[, target] = x
    b.pred = d2 %*% object$coefficients
    p.pred = 1 / (1 + exp(-b.pred))
    drop(p.pred)
}
P.predi(res, "age")

いつものように，代入記号に = を使ったのは，このブログでのバグを避けるためなのでご勘弁を

return は，「文」ではなく「関数」だ!!

お恥ずかしいことながら，私自身，前に一度はまったことを，またやっちゃったので，3 度目はないようにと，自戒の意味を込めて。

関数で，途中で関数を抜ける，しかも特に何の戻り値も持たずにという場合には，引数を持たない return は使っちゃいけないよということ。

foo <- function(bar, baz) {
    if (bar < 0) {
        return
   } else {
       :
   }
}

これは，ダメ。return は 関数なので，引数を持たないで return なんて書くと，return の定義そのものが戻り値になる。従って，「関数から戻る」なんてことは，期待できない。文法エラーにもならず気づかないので，最大限の注意が必要。

書くならば return() とか return(NULL) とか，とにかく，引数を明示的に書くべし。

まあ，面白っちゃあ面白い。そうでなきゃ，クソだ

> fivenum(iris[1:4])
 以下にエラー `[.data.frame`(x, order(x, na.last = na.last, decreasing = decreasing)) :
  undefined columns selected
> fivenum(as.numeric(iris[1:4]))
 以下にエラー fivenum(as.numeric(iris[1:4])) :
   (list) オブジェクトは 'double' に変換できません
> fivenum(as.matrix(iris[1:4]))
[1] 0.1 1.7 3.2 5.1 7.9

ベクトル演算を効率的に使う

“（※ただしイケメンに限る）” をランダムに挿入するR関数 だけど

3点についてコメント

1. 文字レベルで分解する必要はなくて，文に分解するに分解するときに，split に複数の文字を指定するやり方（オンラインヘルプに書いてある）

2. ある確率で何かをやるときは，rbinom ではなく runif を使う

3. for ではなくベクトル演算を行う

ということで，以下のような単純なプログラムになりました。

tadasi2 = function(text, p) {
    x = unlist(strsplit(text, "[。.！]"))
    y = ifelse(runif(length(x)) , "（※ただしイケメンに限る）", "")
    x = paste(x, y, "。", sep="")
    paste(x, collapse = "")
}
x = "英語なんて言葉なんだ！ こんなものやれば誰だってできるようになる！"
tadasi2(x, 0.3)

プログラムも簡潔に書きたい

二人の相性（％）を瞬時に判定するR関数　だけど，計算を簡単に行うプログラムを書くのが面倒だとちょっと...

このプログラムを入力するのが面倒くさい。

というわけで，以下のように書いてみる。

uranai = function(NAME1, NAME2, PRINT = TRUE) {
    ### 名前を統合 ###　
    x = unlist(strsplit(c(NAME1, NAME2), split = ""))
    y = numeric(length(x))
    y[x %in% unlist(strsplit("あぁかさたなはまらがざだばぱやゃわ", ""))] = 1
    y[x %in% unlist(strsplit("いぃきしちにひみりぎじぢびぴ",      ""))] = 2
    y[x %in% unlist(strsplit("うぅくすつぬふむるぐずづぶぷゆゅっ", ""))] = 3
    y[x %in% unlist(strsplit("えぇけせてねへめれげぜでべぺ",      ""))] = 4
    y[x %in% unlist(strsplit("おぉこそとのほもろごぞどぼぽよょ",   ""))] = 5
    ###　計算 ###　
    for (i in 3:length(x)) {
        if (PRINT) {
            print(y)
        }
        n = length(y)
        if (n == 2) {
            break
        }
        y = (y[-n] + y[-1])%%10
    }
    ###　パーセンテージで出力 ###
    cat(paste("ふたりの相性は、", 10 * y[1] + y[2], "%", " です。", sep = ""))
}
uranai("やべひろゆき", "あおきゆうこ")
uranai("うちだゆうや", "きききりん", PRINT = FALSE)
uranai("たむらゆう", "ふくたじゅんや", PRINT = TRUE)

ところで，元のプログラムで

x <- X[i]
ifelse( x == "あ" | x == "か"| ... | x == "ば" ,y[i] = 1,
ifelse(x == "い"...)...)

というのは，

y[i] = ifelse( x == "あ" | x == "か"| ... | x == "ば" , 1,
ifelse(x == "い"...)...)

などとすべし。

まあ，とにかく，ハッピー・ハッキング！！

 

ダミー変数の標準化偏回帰係数

重回帰分析で標準化回帰係数βを出力するR関数　が紹介されている

ちょっと，見にくいところと余分なことがあるので，修正してみた（代入記号に = を使うのは，不本意だが，許して欲しい）

lm.Beta2 = function(res) {
    ysd = sd(res$model[, 1]) # 結果変数のSD
    idv = res$model[, -1, drop = FALSE]
    N = ncol(idv)
    res.beta = sd = res$coefficients[-1]
    for (j in 1:N) {
        xxx = idv[, j]
        if (class(xxx) == "factor") {
            for (i in 2:nlevels(xxx)) {
                lab = paste(colnames(idv)[j], levels(xxx)[i], sep = "") # 対象の変数
                dummy = as.integer(xxx) == i # 1/0 の変数化
                sd[lab] = sd(dummy)
            }

        } else { # 数値だったら簡単にβを計算できる
            lab = colnames(idv)[j] # 対象の変数
            sd[lab] = sd(xxx)
        }
    }
    res.beta * sd / ysd
}

R の他の関数のように，下請けの関数を組み合わせて目的を達成するという方針でやるならば，以下のように簡単に書くことができる（1行にすることもできるが，そこまでは...）

lm 関数が返すオブジェクト object に対して，

lm.Beta3 = function(object) {
  d = model.matrix(object$terms, eval(object$model, parent.frame()))
  object$coefficients[-1] * apply(d[, -1, drop=FALSE], 2, sd) / sd(object$model[,1])
}

> lm.Beta3(res)
        x        zB        zC
0.5238095 0.2436580 0.5685352

楽しい，ハッキングだった(^_^)

ありがとうございました

必要なことと余分なことの区別

上三角行列のデータから対象行列を作る

ちょっと余分なことをしているので，

dist.mat2 = function(x, n) {
    m = matrix(0, n, n)
    m[lower.tri(m)] = x
    return(m+t(m))
}
dist.mat2(1:6, 4)

このブログで代入記号を使うとまずいことになるので = を使う