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),それぞれの長所短所があるのでしょうかね。
# 直接コメントを書こうとしたのですが,どうやったらよいか見つけることが出来ませんでしたので。
言語によって,gsub の仕様が少し変わるが基本は同じ。
BEGIN {
str = "1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989"
gsub("(..)", "&,", str)
print gsub("00|01|04|15|17|18|24", "@", str)
}
(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))))
計算時間を贅沢に
> 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
:
文字列操作で書いてみる(題意通りの解が得られる区切り方)
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年猶予
論文訂正に1年の猶予というのではなく,まずは取り消し,再提出されたならば厳正に審査するというのが正しい対処法ではないか?
なぜ猶予刑なのか。理解に苦しむ。
> 10/06 朝の体重78.6㎏、体脂肪率31%でした。
> 10/05 朝風呂後の体重78.6㎏、体脂肪率27%でした。
> 10/04 朝の体重79.4㎏、体脂肪率26%でした。
体の中の脂肪の量って,そんなに変動するものなのか?
体脂肪量を表示する体重計の精度とか理論とか考えるべきでは?
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 よりは短いしトリッキーでもない