CodeIQ

https://codeiq.jp/challenge.php?challenge_id=893

【問題】
1からnまで連続する正の整数があります。それらの中に、「7」がいくつあるかを数えてください。

【例】
n = 99とすると、
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99

「7」を含む値は、
7 17 27 37 47 57 67 70 71 72 73 74 75 76 77 78 79 87 97

の19個ですが、77には「7」が2つ含まれています。ですので、「7」の数は19個ではなく、20個と答えてください。

n=99
n=77777
n=23678947
n=732465890
n=1912478368
のときの個数を求めてください

締め切り　2014年7月7日 AM8:00

ネタばらし禁止だそうで，締め切り後に，私のプログラムを提示します。

無駄をなくし，短くすると，速くなる（本当か？）

「コーパス間で頻度差の大きい単語を特定する」だが...
http://t.co/aTOkCabdnT

少なくとも，mat3 への代入文は for ループの外に出さないといけない
（ソースフォーマットが悪い（ない）ので，自分で気づけない）
fisher.bat = function(x) {
    sum = apply(x, 2, sum)
    nr = nrow(x)
    nc = ncol(x)
    mat = matrix(0, nr, 1)
    mat2 = cbind(x, p.value = mat)
    for (i in 1:nr) {
        temp = rbind(x[i, ], sum - x[i, ])
        mat2[i, nc + 1] = fisher.test(temp)$p.value
        mat3 = mat2[sort.list(mat2[, nc + 1]), ]
    }
    list(results = mat3)
}

for をなくすと共に，いろいろと無駄をなくしてみる
fisher.bat2 = function(x) {
    sum = colSums(x)
    y = t(sum - t(x))
    p.values = apply(cbind(x, y), 1, function(temp) fisher.test(matrix(temp, 3))$p.value)
    cbind(x, p.values)[order(p.values), ]
}

速くなったようだ
> system.time(fisher.bat(nouns))
   ユーザ   システム       経過  
     0.449      0.017      0.463
> system.time(fisher.bat2(nouns))
   ユーザ   システム       経過  
     0.165      0.012      0.175

ブルートフォースだけじゃダメよ

この下の問題だけど，よく解析してからプログラムを書くと，あっと驚く速さになる。

まずは，探索範囲を狭めること。

> system.time({
+     euclid = function(m, n) {
+         while ((temp = n%%m) != 0) {
+             n = m
+             m = temp
+         }
+         return(m)
+     }
+     N = 1000
+     mx = floor(sqrt((1 + sqrt(2)) * N/2))
+     d = expand.grid(1:mx, 1:mx)
+     d = d[d[, 1] > d[, 2], ]
+     n = d[, 1]
+     m = d[, 2]
+     a = n^2 - m^2
+     b = 2 * n * m
+     c = n^2 + m^2
+     d = cbind(ifelse(a < b, a, b), ifelse(a > b, a, b), c)
+     d = d[d[, 2] < N, ]
+     GCD = apply(d, 1, function(x) Reduce(euclid, x))
+     d = d/GCD
+     d = unique(d)
+ })
   ユーザ   システム       経過  
     0.027      0.002      0.036

ベクトル演算は強力だなぁ

「ピタゴラスの定理の整数解（原始ピタゴラス数）」ですが...

元のプログラムは，私のマシンだと1秒ほどかかる。
> system.time({
+     ## ピタゴラスの定理
+     library(gmp)
+     n = 1000
+     i = 1
+     ans = matrix(0, floor(n) * 3, ncol = 3)
+
+     for (a in 1:n) {
+         for (b in 1:n) {
+             if (a <= b) {
+                 if (floor(sqrt(a * a + b * b)) == sqrt(a * a + b * b)) {
+                     temp = c(a, b, sqrt(a^2 + b^2))
+                     if (Reduce(gcd, temp) == 1) {
+                         ans[i, ] = c(a, b, sqrt(a^2 + b^2))
+                         i = i + 1
+                     }
+                 }
+             }
+         }
+     }
+
+     ans = ans[apply(ans, 1, sum) != 0, ]
+ })
   ユーザ   システム       経過  
     1.004      0.007      1.003

あまりよいプログラムではないが，ベクトル演算を使うと3.7倍ほど速くなる。
> system.time({
+     library(gmp)
+     n = 1000
+     d = expand.grid(1:n, 1:n)
+     s = sqrt(rowSums(d^2))
+     d = cbind(d, s)
+     d = d[floor(s) == s & d[, 1] <= d[, 2], ]
+     GCD = apply(d, 1, function(x) Reduce(gcd, x))
+     d = unique(d/GCD)
+     d = d[order(d[, 1]), ]
+ })
   ユーザ   システム       経過  
     0.280      0.010      0.301

update の価値が理解できないなあ

「複数の学習データから行う回帰分析」なんだけど...

コメントの押し売りをするつもりはないんだけど，コメントの受付方法法くらい用意しておく方がよいのではないかと思います。

筆者は，update という関数名からなのか，「update するたびに素晴らしい結果が約束される」というような幻想を持っているのかな？

update は，以前に行った解析の主として formula を受け継いで，再度分析するということでは？再度の分析が必要なのは，データがアップデートされた（間違いが修正された，観察が増えた）ときに，最新のデータで分析結果をアップデートするということではないか？

例えば，以下のような例。

その時点で得られているデータに基づいて，回帰分析を行う。

> (d <- data.frame(x=c(2,1,3,4,5), y=c(3,5,4,7,9)))
  x y
1 2 3
2 1 5
3 3 4
4 4 7
5 5 9
> lm1 <- lm(y ~ x, data=d)
> summary(lm1)

Call:
lm(formula = y ~ x, data = d)

Residuals:
   1    2    3    4    5
-1.4  1.8 -1.6  0.2  1.0

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   2.0000     1.7963   1.113    0.347
x             1.2000     0.5416   2.216    0.114

Residual standard error: 1.713 on 3 degrees of freedom
Multiple R-squared:  0.6207,    Adjusted R-squared:  0.4943
F-statistic: 4.909 on 1 and 3 DF,  p-value: 0.1135

新たにデータが加わった。（同じ名前でもよいが，データフレームの名前を付け替えた）

> (d2 <- rbind(d, c(6, 12)))
  x  y
1 2  3
2 1  5
3 3  4
4 4  7
5 5  9
6 6 12

新たなデータを使って，formula はそのまま（lm1 に記録されているので，そこから取り出して使う）

> summary(update(lm1, data=d2))
Call:
lm(formula = y ~ x, data = d2)

Residuals:
       1        2        3        4        5        6
-1.26667  2.33333 -1.86667 -0.46667 -0.06667  1.33333

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   1.0667     1.6479   0.647   0.5527
x             1.6000     0.4231   3.781   0.0194

Residual standard error: 1.77 on 4 degrees of freedom
Multiple R-squared:  0.7814,    Adjusted R-squared:  0.7267
F-statistic:  14.3 on 1 and 4 DF,  p-value: 0.01941

しかし，その結果は，最新のデータフレームを使って，lm1 を得たときと同じ formula を使って lm 関数により解を求めたものと，全く同じなのだ。

> summary(lm(y ~ x, data=d2))

Call:
lm(formula = y ~ x, data = d2)

Residuals:
       1        2        3        4        5        6
-1.26667  2.33333 -1.86667 -0.46667 -0.06667  1.33333

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   1.0667     1.6479   0.647   0.5527
x             1.6000     0.4231   3.781   0.0194

Residual standard error: 1.77 on 4 degrees of freedom
Multiple R-squared:  0.7814,    Adjusted R-squared:  0.7267
F-statistic:  14.3 on 1 and 4 DF,  p-value: 0.01941

formula が長ったらしくたって，それを記録しておいて再度使えば良いだけだし。

lm オブジェクトとして保管しておき，その中の formula を使うなどということは，あまりメリットはないように思うが？？？？

余計なことはしない

「統計を始める方へ(1)_データ環境Rの基本的なプログラミング｜データアーティスト」で，ついでに...

z <- c(1:n) のようなプログラム片が示されていたけど，c 関数は余計

> n <- 1e8
> system.time(z <- c(1:n))
   ユーザ   システム       経過  
     0.431      0.297      0.724
> system.time(z <- 1:n)
   ユーザ   システム       経過  
     0.112      0.152      0.261

最適な関数を最適に使おう

「統計を始める方へ(1)_データ環境Rの基本的なプログラミング｜データアーティスト」で，サイコロをシミュレーションするプログラム片で，round(6*runif(n)+0.5) が使われていたので，重箱の隅つつき...

> n <- 1e8
> system.time(x <- round(6*runif(n)+0.5))
   ユーザ   システム       経過  
     7.789      0.625      8.362
> system.time(x <- sample(1:6, n, replace=TRUE))
   ユーザ   システム       経過  
     2.860      0.460      3.296
> system.time(x <- sample(6, n, replace=TRUE))
   ユーザ   システム       経過  
     1.780      0.153      1.923

意外とやりそうなのが，sample(1:6, n, replace=TRUE) であるが，これは，sample(6, n, replace=TRUE) より劣る。以下の例を見れば分かる。何故そうなるかは，sample のソースを見れば分かる。

> n <- 1e8
> m <- 1e5
> system.time(y <- sample(1:m, n, replace=TRUE))
   ユーザ   システム       経過  
     3.281      0.485      3.741
> system.time(y <- sample(m, n, replace=TRUE))
   ユーザ   システム       経過  
     1.781      0.162      1.932