ベンフォードの法則の例

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)

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