ベンフォードの法則の例

2のべき乗数，2^0, 2^1, 2^2, 2^3, ... の先頭桁を考える。

一見，1, 2, 4, 8, 1, 3, 6, 1, 2, 5, 1, 2, 4, ... と長さ10の周期で循環しそうに見えるがそれは違う。

先頭桁の分布はベンフォードの法則に従う。

> library(gmp)
> func <- function(n)
+ {
+     a <- as.character(n)
+     a <- strsplit(a, "")
+     a <- unlist(a)
+     return(a[1])
+ }
> a <- integer(100000)
> for (i in 1:100000) {
+     n <- as.bigz(2)^i
+     a[i] <- func(n)
+ }
> table(a)

このプログラムによって得られる結果は，正確にベンフォードの法則に従うものである。

出現度数は次のようになる。

a
    1     2     3     4     5     6     7     8     9
30102 17611 12492  9692  7919  6695  5797  5116  4576

3 のべき乗の場合も同じような結果になる（驚くほど似ている）。

a
    1     2     3     4     5     6     7     8     9
30101 17611 12492  9692  7917  6696  5799  5116  4576

いずれの場合も，各数字の理論的な出現確率は

[1] 0.30103000 0.17609126 0.12493874 0.09691001 0.07918125 0.06694679
[7] 0.05799195 0.05115252 0.04575749

なのだ。

 

みったくなすの ggplot2

散布図も

p <- ggplot(iris, aes(Sepal.Width, Sepal.Length))
p2 <- p + geom_point(aes(colour = Species))
p3 <- p2 + xlab("Sepal width") + ylab("Sepal length") +
            xlim(2, 4.5) + ylim(4, 8)
print(p3)

old <- par(mar=c(3, 3, 1, 1), mgp=c(1.8, 0.6, 0), las=1)
plot(iris$Sepal.Width, iris$Sepal.Length, col=c(1, 2, 4)[as.integer(iris$Species)], pch=19, cex=0.7, xlab="Sepal width", ylab="Sepal length")
legend("topright", legend=c("setosa", "versicolor", "virginica"), col=c(1, 2, 4), pch=19, cex=0.8, bty="n")
par(old)

どっちがきれいか，明らかだと思うが

 

ggplot2 をもてはやす輩がいるようだが

たいして，きれいでもないし，大局観に基づいているわけでもないし，ブラックボックス信奉を煽るだけではないかなあ。そのうち，下火になるだろう。

背景がグレー（デフォルトがグレーなだけで，背景を白にすることもできるけど）というのは，いかにもださい。色も「落ち着いた色」を使おうとしているのだが，くすんだ色で，ださい感じ。鮮やかな色というわけではないけど，いかにもみすぼらしい色使いだ。

格子線を描くのもださい。まあ，そのあたりの感覚は個人差があるのだろうけど。まあ，Excel のださいグラフに通じるところがあるという気がするわけだ。あれも，描画領域がグレーで，補助線引きマクリだったよね。そういうのが好きな人は好きなんだろうね（ちなみに，表を作るときにセルの周りにガチガチに線を引き回った表を作る習慣ができたのははエクセルのせいだよね。古き良き時代の書籍における表はキレイだったけど，最近の本にはエクセルで作りましたというガチガチに線を引いた表がはびこってるね）。

フォントがださい。Windows 環境なら，もっとださい。日本語描けないなんて，英語が不得意な私にとっては，一番のマイナス点だね。それだけで，候補から堕ちる。

レイヤーを重ねるなんていっても，所詮は，描画（効果）を累積するだけであって，それを対話的にやるなんてことは現実的ではないわけで，プログラムする訳なんだから，結局はできあがったものだけを評価するんだから，途中を + で繋ごうが，別々の関数を呼ぼうが大差はない。

同じ効果を持つ図を描こうとすれば，それなりの関数を定義すればよいだけだ（何十種類ものグラフを描く必要のある人はいないだろう）。

前半の図

library(ggplot2)
qplot(wt, mpg, data = mtcars) +
  stat_smooth(method = "lm", se = FALSE)

後半の図

old <- par(mar=c(2.5, 2.5, 1, 1), mgp=c(1.6, 0.6, 0))
plot(mtcars[c("wt", "mpg")], pch=19, cex=0.7, las=1)
abline(lm(mpg~wt, mtcars), col="blue")
par(old)

ブランド志向はここまでくるか

> 依然としてよくわからないので間瀬先生にメールにて質問させていただいたところ、即座に親切なご回答をいただきましたので紹介しますと、”「正確に」答え られる概念ではありませんし、メイルで答えられるような簡単な話でもありません。マルコフ連鎖モンテカルロ（MCMC）法に関する文献を探してください。 詳しい議論があるでしょう。”

「正確に」答え られる概念ではありませんし、メイルで答えられるような簡単な話でもありません。マルコフ連鎖モンテカルロ（MCMC）法に関する文献を探してください。 詳しい議論があるでしょう。

というのは，それ以前に得られているコメントと内容的には同じだよね。

まあ，違う所をしいてあげれば，「質問者のプライドを傷つけない」というところかな。

質問者も，ナイーブになったものだ（念のため，ナイーブというのは，英語圏では負のイメージを持つのだよ）

 

 

n日後を返す関数を返す関数

http://ja.doukaku.org/29/

整数nを渡すと「日時のデータを受け取って、n日後の日時を返す関数」を返す関数を作ってください。関数を返す関数が作れない場合は、関数の代わりになるようなオブジェクトでも構いません。

Pythonの対話的インタプリタで表現すると下のようになります。

>>> def n_days_later(n):
        ?????

>>> five_days_later = n_days_later(5)
>>> datetime.datetime.now()
datetime.datetime(2007, 7, 20, 20, 11, 42, 78000)
>>> five_days_later(_)
datetime.datetime(2007, 7, 25, 20, 11, 42, 78000)

出題の意図としては「関数を返す関数」と「日時の差分の扱い方」を、それぞれ単体だと簡単すぎるので合わせ技にしてみた、というところです。

「2038年問題に注意せよ」と言う注意書きがあるけど

func <- function(s=Sys.Date(), n=0)
{
    substr(as.POSIXct(s)+n*86400, 1, 10)
}

実行例

> func()
[1] "2011-02-17"
> func("2035-01-01", 10000)
[1] "2062-05-19"

与えられた数字のケタ数

http://ja.doukaku.org/40/

与えられた数字のケタ数と、最大桁の位を求めてください。
数字が2469なら4桁で最大桁は1000の位です。
600なら3と100、1なら1と1です。

言語によっては，簡単すぎる問題もある。

func <- function(n)
{
    k <- length(unlist(strsplit(as.character(n), "")))
    return(c(k, 10^(k-1)))
}

実行例
> func(2469)
[1]    4 1000
> func(600)
[1]   3 100
> func(1)
[1] 1 1

bigq クラスの数学関数

sqrt はまあ実用的

sqrt.bigq <- function(x, eps=1e-50)
{
    x2 <- as.bigq(sqrt(as.double(x)))
    repeat {
        x3 <- (x2+x/x2)/2
        if (abs(x2-x3) < eps) break
        x2 <- x3
    }
    return(x3)
}
# as.double(sqrt.bigq(5000.45789))^2

# sqrt(as.bigq(3.456))

3 乗根

qubic.root <- function(x, eps=1e-50)
{
    x2 <- as.bigq(as.double(x)^(1/3))
    repeat {
        x3 <- (x/(x2*x2)+2*x2)/3
        if (abs(x2-x3) < eps) break
        x2 <- x3
    }
    return(x3)
}
# as.double(qubic.root(123.4566789))^3

exp は取りあえず使えればよいものということで

exp.bigq <- function(x, eps=1e-50, n.term=1000)
{
    ans <- xp <- fac <- as.bigq(1)
    for (i in 1:n.term) {
        xp <- xp*x
        fac <- fac*i
        old <- ans
        ans <- ans+xp/fac
        if (abs(old-ans) < eps) {
            return(ans)
        }
    }
    warnings("収束せず")
    return(ans)
}
# log(as.double(exp.bigq(40.235)))

# exp(as.bigq(1.234))

 

 

文字列型日時ののN秒後時間取得

http://ja.doukaku.org/204/

日時を表す文字列と時間(秒)を受け取り
指定された日時からN秒となる日時を出力する関数 DateEx() を作成してください。

関数の仕様は次の通りです。
1. 入力となる日時の書式は任意である。
    → プログラムの都合に合わせてよい。
2. 入力となる時間(秒)は、負の値も許容すること。
    また、負の値が指定された場合、指定の日時よりも前の日時を出力すること
3. 出力する日時は入力の日時と同じ書式をとる文字列であること
4. 出力する日時は正規化されていること
5. 出力先は標準出力、または、バッファのいずれでもよい。

たとえば、DateEx("20080827235925",40)ならば
出力は
「20080828000005」です。

余力があれば時間を省略可能とし、
省略された場合は「現在時刻」を利用するようにしてみてください。

DateEx <- function(str=Sys.time(), N=0)
{
	return(as.POSIXct(str)+N)
}

実行例
> DateEx("2008-08-27 23:59:25 JST", 40)
[1] "2008-08-28 00:00:05 JST"
> DateEx() # 現在時刻
[1] "2011-02-16 23:28:43 JST"

回文日付

http://ja.doukaku.org/289/

1. YYYYMMDD形式で書き下した日付が回文になる日を探して列挙してください。探すのは西暦4桁の範囲（YYYY年MM月DD日）とします。

例えば、2010年1月2日をYYYYMMDD形式で書き下すと 20100102 となり、逆から読んでも 20100102 になります。

また、日付を扱うライブラリを利用した場合はそのように明記してください。

余裕のある方は以下にもチャレンジしてください。

1. 他の形式（e.g. DDMMYYYY形式など）で回文になる日付を探せるように改良してください。

※ちなみに、日本では年月日の順で表記するのが一般的ですが、表記の順番は国によってバラバラのようです。

R はベクトル演算が得意

d <- seq(as.Date("1000-01-01"), as.Date("9999-12-31"), by=1)
d1 <- matrix(unlist(strsplit(as.character(d), "")), 10)
d2 <- (d1[1,] == d1[10,]) & (d1[2,] == d1[9,]) & (d1[3,] == d1[7,]) & (d1[4,] == d1[6,])
print(d[d2])

331 個ある。古い White MacBook で 84 秒ほどかかる。

> system.time({
+ d <- seq(as.Date("1000-01-01"), as.Date("9999-12-31"), by=1)
+ d1 <- matrix(unlist(strsplit(as.character(d), "")), 10)
+ d2 <- (d1[1,] == d1[10,]) & (d1[2,] == d1[9,]) & (d1[3,] == d1[7,]) & (d1[4,] == d1[6,])
+ print(d[d2])
+ })
  [1] "1001-10-01" "1010-01-01" "1011-11-01"
  [4] "1020-02-01" "1021-12-01" "1030-03-01"

：

[325] "9230-03-29" "9240-04-29" "9250-05-29"
[328] "9260-06-29" "9270-07-29" "9280-08-29"
[331] "9290-09-29"
   ユーザ   システム       経過 
    84.676     16.192    172.650

バイナリクロック

http://ja.doukaku.org/275/

時刻を二進数相当の表現で出力する時計アプリケーションを書いてください。
20:18の場合，例えば以下の様な出力をするイメージです。

出力例:
■□■□□
□■□□■□

グラフィックウインドウに描くことにしよう

func <- function(a, s, e, y)
{
    text(1, y, paste(ifelse(rev(rawToBits(as.raw(substr(a, s, e)))), "x", "o"), collapse=""))
}
plot(c(0, 2), c(0, 5), type="n", xlab="", ylab="", xaxt="n", yaxt="n", bty="n")
repeat {
    a <- Sys.time()
    rect(0.5, 0.5, 1.5, 4.5, col="white")
    text(1, 4, substr(a, 12, 19))
    func(a, 12, 13, 3)
    func(a, 15, 16, 2)
    func(a, 18, 19, 1)
    repeat {
        if (substr(a, 18, 19) != substr(Sys.time(), 18, 19)) break
    }
}

文字の出現頻度？

まあ，これはこれでご立派なのだが

hash <- new.env(hash=TRUE, parent=emptyenv())
word <- 'abracadabra'
for (c in strsplit(word, NULL)[[1]]) {
  if (exists(c, hash)) {
    assign(c, get(c, hash) + 1, hash)
  } else {
    assign(c, 1, hash)
  }
}

for (c in ls(hash)) {
  cat(sprintf ('%s : %d\n', c, get(c, hash)))
}

word <- 'abracadabra'
table(strsplit(word, NULL)[[1]]) で十分でしょう

> word <- 'abracadabra'
> table(strsplit(word, NULL)[[1]])

a b c d r
5 2 1 1 2

メンバーシップテストって何だ？

vec <- vector()
for (x in unlist(strsplit("abracadabra", NULL))) {
  if (! x %in% vec) {
    vec <- append(vec, x)
  }
}
print (vec)

[1] "a" "b" "r" "c" "d"

重複を除いた要素を得ることか？

unique(unlist(strsplit("abracadabra", NULL))) が一番簡単。

> unique(unlist(strsplit("abracadabra", NULL)))
[1] "a" "b" "r" "c" "d"

これだと，出現単語表を作ったりすることもできる。

> z <- c("one", "one", "can", "can", "too", "bad")
> unique(z)
[1] "one" "can" "too" "bad"

ベクトルの最後の要素

コメントされるのがイヤならば，コメント不許可にしておけばよい。せっかくなのでこちらに書いておく。

> system.time(for (i in 1:100000) {n <- length(x); x[n]})
   ユーザ   システム       経過 
     0.146      0.002      0.147
> system.time(for (i in 1:100000) tail(x, n=1))
   ユーザ   システム       経過 
     4.236      0.020      4.223

tail の中を見てみれば分かるが，いろいろエラー処理などをやっているが，所詮は x[length(x)] をやることになっている。いろいろな余分な処理のためひどく時間が掛かっている。

もう一つの案が次に示すもの。ベクトルを反転させているので時間が掛かる。長いベクトルならそれだけ余計な時間が掛かる。でも，まあ，たいした時間ではない。

> system.time(for (i in 1:100000) rev(x)[1])
   ユーザ   システム       経過 
     1.486      0.007      1.481

全角数字列を半角数字列へ

encoding 依存かな。Mac なら以下も。

> x <- "２００８"
> paste(as.character(as.integer(sapply(unlist(strsplit(x, "")), charToRaw)[3,])-9*16), collapse="")
[1] "2008"

> paste(as.character(as.integer(sapply(unlist(strsplit(x, "")), charToRaw)[3,])%%16), collapse="")
[1] "2008"

 

Windows ならば

paste(as.character(as.integer(sapply(unlist(strsplit(x, "")), charToRaw)[2,])-79), collapse="")

 

データフレームの左右をひっくり返す

> (a <- data.frame(x=1:5, y=letters[1:5], z=rnorm(5)))
   x y                  z
1 1 a -0.06277869
2 2 b  0.99328294
3 3 c -1.00593094
4 4 d -0.75098495
5 5 e  1.74401638
> (b <- rev(a))
                    z y x
1 -0.06277869 a 1
2  0.99328294 b 2
3 -1.00593094 c 3
4 -0.75098495 d 4
5  1.74401638 e 5