効果量

r-statistics-fanの日記
Rで標準化効果量を計算する(Rewriting the code of SAS proc %stddiff to R)
http://r-statistics-fan.hatenablog.com/entry/2014/10/14/215903

比率と，順序尺度については大いに参考になりました。
間隔尺度についての効果量は，各群のサンプルサイズも考慮して
effect size g = (mu1-mu0)/sd
sd = sqrt( (n1-1)sd1^2 + (n0-1)sd0^2) / (n1+n0-2) )
みたいに計算することが多いようですが（Glass 1976），それぞれの長所短所があるのでしょうかね。

# 直接コメントを書こうとしたのですが，どうやったらよいか見つけることが出来ませんでしたので。

ラズベリーは何個ある？（awk）

言語によって，gsub の仕様が少し変わるが基本は同じ。

BEGIN {
str = "1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989"
gsub("(..)", "&,", str)
print gsub("00|01|04|15|17|18|24", "@", str)
}

ラズベリーは何個ある？（3）

（2）と似ているが

str = "1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989"
length(grep("00|01|04|15|17|18|24", unlist(strsplit(gsub("(..)", "\\1,", str), ","))))

また，以下のようにも

nchar(gsub("[^@]", "", gsub("00|01|04|15|17|18|24", "@", gsub("(..)", "\\1,", str))))

2番目以降の四つ子素数

計算時間を贅沢に

> library(gmp)
> n = 30*0:10000
> n = n[isprime(n+11) > 0 & isprime(n+13) > 0 & isprime(n+17) > 0 & isprime(n+19) > 0]
> t(sapply(n, "+", c(11, 13, 17, 19)))
        [,1]   [,2]   [,3]   [,4]
 [1,]     11     13     17     19
 [2,]    101    103    107    109
 [3,]    191    193    197    199
 [4,]    821    823    827    829
 [5,]   1481   1483   1487   1489
 [6,]   1871   1873   1877   1879
 [7,]   2081   2083   2087   2089
 [8,]   3251   3253   3257   3259
 [9,]   3461   3463   3467   3469
[10,]   5651   5653   5657   5659
：

ラズベリーは何個ある？（2）

文字列操作で書いてみる（題意通りの解が得られる区切り方）

str = "1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989"
str = gsub("(..)", "\\1,", str)
nchar(str)-nchar(gsub("00|01|04|15|17|18|24", "*", str))

ラズベリーは何個ある？

> 最初の2桁は 14 なので、アルファベット・インデックス表の 「o」 にあたり、カウントしません。 > 次の2桁は 15 なので、アルファベット・インデックス表の 「p」 にあたり、カウントします。

以下のプログラムは，要求されている仕様と異なる。以下のプログラムでは，14, 41, 15  のようにして二桁ずつ区切る。ははは。

この解は，βバージョン

str = "1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989"
str1 = unlist(strsplit(str, ""))
str2 = str1[2:1000]
str1 = str1[1:999]
func = function(x, y) {
    (x == 0 && y == 0)  || # a
    (x == 0 && y == 1)  || # b
    (x == 0 && y == 4)  || # e
    (x == 1 && y == 5)  || # p
    (x == 1 && y == 7)  || # r
    (x == 1 && y == 8)  || # s
    (x == 2 && y == 4)     # y
}
    
system.time(print(sum(mapply(func, str1, str2))))

[1] 63
   ユーザ   システム       経過  
     0.030      0.001      0.037

題意どおりに解く場合は，str1, str2 をちょっと変えればよい。答えは 31 になるはず。hahaha...

白か黒か，はっきりせよ!

> 早大:小保方氏の博士号取り消し決定　論文訂正に１年猶予

論文訂正に1年の猶予というのではなく，まずは取り消し，再提出されたならば厳正に審査するというのが正しい対処法ではないか？

なぜ猶予刑なのか。理解に苦しむ。

何が測定されているのか分かっているのだろうか

> 10/06 朝の体重78.6㎏、体脂肪率31％でした。
> 10/05 朝風呂後の体重78.6㎏、体脂肪率27%でした。
> 10/04 朝の体重79.4㎏、体脂肪率26％でした。

体の中の脂肪の量って，そんなに変動するものなのか？

体脂肪量を表示する体重計の精度とか理論とか考えるべきでは？

Pumpkin Helmet

R の次のバージョンのコードネームは Pumpkin Helmet だそうだ。

いつものように，Snoopy と Charlie Brown ゆかりのもの。

http://dvd.netflix.com/Movie/You-re-a-Good-Sport-Charlie-Brown/70111491

最初の絵の説明中には，以前のコードネーム Good Sport(R 3.0.1) や Masked Marvel(R 3.0.0) も出ている。

http://en.wikipedia.org/wiki/You%27re_a_Good_Sport,_Charlie_Brown も参照。

条件を満たす素数の組

outer は，実に使いでのある関数だ

# 差が 6 の素数の組
n = 2:500
m = n[1-outer(n, n)]
index = which(outer(m,  m, "-") == 6, arr.ind=TRUE)
cat(paste(apply(matrix(m[index], ncol=2), 1, function(x) sprintf("(%i,%i)", x[2], x[1])), collapse=","))
# または
cat(paste(mapply(function(x, y) sprintf("(%i,%i)", x, y), m[index[,2]], m[index[,1]]), collapse=","))
出力
(5,11),(7,13),(11,17),(13,19),(17,23),(23,29),(31,37),(37,43),(41,47),(47,53),(53,59),(61,67),(67,73),(73,79),(83,89),(97,103),(101,107),(103,109),(107,113),(131,137),(151,157),(157,163),(167,173),(173,179),(191,197),(193,199),(223,229),(227,233),(233,239),(251,257),(257,263),(263,269),(271,277),(277,283),(307,313),(311,317),(331,337),(347,353),(353,359),(367,373),(373,379),(383,389),(433,439),(443,449),(457,463),(461,467)

数列の和

1:N までの数列から，2の倍数(2, 4, 6, ...)と3の倍数(3, 6, 9, ...)を取り除いた数列の和

func1 = function(N) {
    n = 1000000
    x = 1:n
    x[1:(n%/%2)*2] = 0
    x[1:(n%/%3)*3] = 0
    x = x[x != 0]
    if (N > length(x)) return(NA)
    else return(sum(x[1:N]))
}
func2 = function(N) {
    m = N%/%2
    s = 1
    for (i in seq_len(m)) {
        s = s+12*i
        s = s-(i == m && N%%2==0)*(6*i-1)
    }
    s
}
func3 = function(N) {
    m = N%/%2
    1 + sum(12 * seq_len(m)) - (!N%%2) * (m * 6 + 1)
}

> func3(1)
[1] 1
> func3(2)
[1] 6
> func3(3)
[1] 13
> func3(4)
[1] 24
> func3(5)
[1] 37

> N = 12345
> func1(N)
[1] 228598537
> func2(N)
[1] 228598537
> func3(N)
[1] 228598537

> library(rbenchmark)
> benchmark(func1(N), func2(N), func3(N))
      test replications elapsed relative user.self sys.self user.child sys.child
1 func1(N)          100  10.722   1787.0    10.275    0.624          0         0
2 func2(N)          100   1.383    230.5     1.450    0.013          0         0
3 func3(N)          100   0.006      1.0     0.005    0.001          0         0

群数列の問題だな

func2 <- function(n) {
    stopifnot(n <= 1388888888888885)
    # 桁数の同じ数を一群とする (5,6,7,8,9), (10,11,...,99), (100,101,...,999), (1000,1001,...,9999),...
    head <- 5 # 群の最初の数値 5, 10, 100, 1000,...
    kosuu <- 5 # 群に含まれる数値の個数 5, 90, 900, 9000,...
    keta <- 1  # 群に含まれる数値の桁数 1, 2, 3, 4,...
    end <- 0 # 群の最後の数値の1の位までの文字数（最初を1と数えて）
    repeat {
        begin <- end + 1 # 群の最初の数値の最上位までの文字数 1, 6, 186, 2886,...
        end <- begin + keta * kosuu - 1 # 群の最後の数値の1の位までの文字数（最初を1と数えて）5, 185, 2885, 38885,...
        if (begin <= n && n <= end) { # n がその群に含まれる数値中にある n = 215 なら，100 で始まる 3 桁の数字の群の中にある186 < 215 < 2885
#            dump(c(head=head, kosuu=kosuu, keta=keta, begin=begin, end=end))
            mojiiti <- n - begin # 群の最初の数値の最上位から何番目の文字か 215-186=29番目。ただし，0から数える。
            number <- head + mojiiti %/% keta # どの数値中の文字か 100 + 29 %/% 3 = 109 の中の 29%%3=2番目。ただし，0から数える。
            ans <- unlist(strsplit(as.character(number), ""))[mojiiti %% keta + 1] # 数値中のどの文字か
            return(as.integer(ans))
        }
        head <- 10 ^ keta # 次の群のために更新
        kosuu <- 9 * head
        keta <- keta + 1
    }
}
# head            kosuu         keta    begin             end
#  100            900                3    186                 2885
#  1000            9000            4    2886             38885
#  10000        90000            5    38886             488885
#  100000        900000            6    488886             5888885
#  1000000        9000000            7    5888886             68888885
#  10000000        90000000        8    68888886         788888885
#  100000000        900000000        9    788888886         8888888885
#  1000000000    9000000000        10  8888888886         98888888885
#  10000000000        90000000000        11    98888888886        1088888888885
#  100000000000        900000000000        12    1088888888886    11888888888885
#  1000000000000    9000000000000    13  11888888888886    128888888888885
#  10000000000000    90000000000000    14  128888888888886 1388888888888885
options(scipen=10)
func2(579334) # 7
func2(62346205) # 6
func2(787925698) # 9
func2(1234567890) # 5
func2(3333333333) # 0
func2(8888888886) # 1
func2(58888886) # 4
func2(186) # 1
func2(2886) # 1
func2(38886) # 1
func2(488886) # 1
func2(5888886) # 1
func2(68888886) # 1
func2(788888886) # 1
func2(9888888886) # 1
func2(1088888888886) # 1
func2(11888888888886) # 1
func2(1388888888888885) # 9

多倍長足し算（引き算も）

負の数は10の補数とするので，引き算もできる

N <- 31 # 符号を含めた桁数
Add <- function(a, b) {
    carry <- function(s) {
        for (i in 1:N) {
            if (s[i] > 9) {
                s[i+1] <- s[i+1]+1
                s[i] <- s[i]-10
            }
        }
        return(s)
    }
    decomp <- function(s) {
        if (sign <- grepl("^-", s)) s <- sub("^-", "", s)
        s <- rev(tail(c(rep(0, 31), as.integer(unlist(strsplit(s, "")))), 31))
        if (sign) {
            s <- 9-s
            s[1] <- s[1]+1
            s <- carry(s)
        }
        return(s)
    }
    c <- carry(decomp(a)+decomp(b))
    sign <- ""
    if (c[N] == 9) {
        sign <- "-"
        c <- 9-c
        c[1] <- c[1]+1
        c <- carry(c)
    }
    length(c) <- N
    c <- paste(sign, sub("^0+", "", paste(rev(c), collapse="")), sep="")
    if (c == "") c <- "0"
    return(c)    
}

> Add("1234", "3456")
[1] "4690"
> Add("1234", "-3456")
[1] "-2222"
> Add("-1234", "3456")
[1] "2222"
> Add("-1234", "-3456")
[1] "-4690"
> Add("99999999", "99999997")
[1] "199999996"

パスカルの三角形

R だと

invisible(sapply(1:10, function(n) cat(choose(n, 0:n), "\n")))

for を使う方が短いが

for (i in 1:10) cat(choose(i, 0:i), "\n")

choose なんていう関数を使うなということならば

x <- 1; for (i in 1:10) cat(x <- c(x, 0)+c(0, x), "\n")

いずれにしろ，Ruby よりは短いしトリッキーでもない