裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

なんで,こんな変な物を好むのだろうか?

2015年01月30日 | ブログラミング

条件に合う行数をカウントする方法」なんだけど,いろいろやってもうまくいかないことがあるのでと愚痴った挙げ句の果てに,

> library(dplyr)
> count(iris, Petal.Length>5)

Source: local data frame [2 x 2]
   Petal.Length > 5 n
1             FALSE 108
2              TRUE 42

が,「dplyr は直感的だし,とても便利」なんですと。私の Mac じゃ,動かなかったけど。

> table(iris$Petal.Length>5)

FALSE  TRUE
  108    42

で十分じゃん!!基本機能でできるんだから。

足らぬ足らぬは,工夫が足らぬ(いや,基本的な知識が足らぬ)

> lengthを駆使して、length(XXX[XXX$xxx>0,])としても何故か欲しい値が返ってきません。

そんなことしても,欲しい値が帰ってこないのはあたりまえ。nrow を使いなさいな。

> nrow(iris[iris$Petal.Length>5,])
[1] 42

または,もっと簡単に(というか,そうするのが普通だが),以下のように。

> sum(iris$Petal.Length>5)
[1] 42

何事も中途半端なら,ほかの何をやっても中途半端。

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

経路は何通り?(その2)

2015年01月29日 | ブログラミング

ブルート・フォースで歯が立たないものは,ちゃんと考えれば,簡単に答が出る。

Excel なんかを立ち上げて,再帰的にシコシコ計算式を書き込んでいけば,数分で答が出る。

R に直すのも,簡単。解答例は,この記事のコメントを参照。

n = 20 のとき,所要時間 0.1 秒以下で,3534526380 を得る(2 ≦ n ≦ 19 のときの解も,その他のノードの解も同時に求まってしまう)。

コメント (2)
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

経路は何通り?

2015年01月28日 | ブログラミング

ローマ水道は、縦横に流れる水が合流する地点のいくつかにおいて橋をかけ、立体交差させることで、水の流れが衝突しないように工夫されています。
下の図は、ローマ水道をシンプルにしたものです。


縦と横にn区間ずつある格子状に水道があるとして、A地点からB地点に最短距離で移動したいのですが、左下と右上を結ぶ対角線上の点のうち、スタートとゴール以外の水の流れが衝突する箇所には橋をかけ、立体交差させています。
つまり、この交差点では曲がることができません。
※図では立体交差を表現するために道を曲げていますが、立体交差の上を超えても下をくぐっても距離は変わらないものとします。

図のようにn=4のとき、スタートからゴールまでの道順は28通りあります。
では、n=10のときとn=20のときについて、スタートからゴールまでの道順が何通りあるかを求めてください。

n = 10 のときは無理矢理求めた(この記事のコメントを参照)が,そのやり方は n = 20 では通用しない。

数列の法則も見いだせない(なさそうに思う)。

コメント (2)
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

とっても大きな数の下 6 桁

2015年01月28日 | ブログラミング
3F1000000 (3 の フィボナッチ数列の第100万項乗)の下 6 桁を求めて下さい。念のため,
3F10 = 355 = 174449211009120179071170507 この数の下6桁は 170507

前にも,同様の質問があったけど,今回はそのときよりも大きな数になる。
しかし,ひるむ必要はない。

解答例は,この記事のコメントを参照。
コメント (1)
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

コメントをもらえるようにしておくと吉

2015年01月27日 | ブログラミング

Web ページを見ていると,他人事ながら,一言お伝えしたいという気になることがある(そういうのを,余計なお節介というんだけど)。

R で全角数字を半角数字に書き換える」において...

> 郵便局のサイトにある「事業所の個別郵便番号データダウンロード」を

> 配布されているのがzip形式なので、URLで直接読み込むのではなくPC上で解凍してから読み込みます。

> > zip <- read.csv("JADD1407.CSV", header=FALSE)  # ダウンロードしたファイルの読み込み
Error in type.convert(data[[i]], as.is = as.is[i], dec = dec, na.strings = character(0L)) :
  invalid multibyte string at...
> 読み込めない・・・

> 郵便局のサイトの説明によるとShift-JISらしいので、文字コードを指定してみる。「encoding="SHIFT-JIS"」を追加します。
> なぜかエラーで読み込めない・・・

> よく分からないので、テキストエディタでcsvをUTF-8に変換してから読み込むこととします。

追加すべき引数は,fileEncoding そして指定するのは SHIFT-JIS ではなく cp932
つまり,

zip <- read.csv("JADD1407.CSV", header=FALSE, fileEncoding="cp932")

> さてここから、冒頭のページに載ってたgsub()というコマンドで、全角数字を半角数字に置き換えてみます。
> しかしgsub()は、1種類の文字列を1種類の文字列に置き換えるものなので、「0~9」→「0~9」をまとめて指定することができないっぽい。なので、「0」→「0」から順番に繰り返していく関数を作ります。
> また、データフレームに対してgsub()による置換を実行すると謎な出力がでてきてむかついたので、列ごとにベクトルとして処理していくことにします。

> gsub.multi <- function (pattern, replacement, data) {
+    # patternには、置き換える対象となる文字列をベクトルで与える
+    # replacementには、上記のそれぞれを何で置き換えるかをベクトルで与える
+    # dataにはデータフレームを入れる
+   
+    for (i in 1:ncol(data)) {
+       for (j in 1:length(pattern)) {
+       data[,i] <- gsub(pattern[j], replacement[j], data[,i])
+       }
+    }
+    return(data)
+ }
> zenkaku <- c("0","1","2","3","4","5","6","7","8","9")
> hankaku <- c("0","1","2","3","4","5","6","7","8","9")
>
> address2 <- gsub.multi(pattern=zenkaku, replacement=hankaku, data=address)

これでもよいけど,chartr という関数が便利。

> chartr("0123456789", "0123456789", "3450246813579")
[1] "3450246813579"

コメント (2)
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

パネルの塗り分け方法

2015年01月26日 | ブログラミング

R, G, B の 3 文字から重複を許して 5 文字を取り出し,1 列に並べる。隣り合う文字が同じにならない文字列を列挙し,最終行に何通りあるかを出力せよ。

3 文字から重複を許して 3 文字を取り出して...の場合は以下のような出力となるように。というのが仕様

RGR
BGR
RBR
GBR
GRG
BRG
RBG
GBG
GRB
BRB
RGB
BGB
12

解答例はこの記事のコメントを参照

コメント (1)
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

ヒントがベタな暗号を解けと...

2015年01月25日 | ブログラミング

キーワード: HAL コンピュータ

以下の暗号を平文に戻せ

"LZH SN STJHZTMNMZQZ JNMN ATMRGNT VN OTQNFTQZLT CD JZHCNJT RHMZRZH RNQDJTQZH CDJHMZJXZ VZSZRH GZ ETSZQH FZ STJHZT MNVN LHSNLDMZHXN"

mrkBDB(sxD3y4Dpi(gf:ja),sxD3y4Dpi(m(gg:ja,gf)),k) # 49 文字解(別の暗号化)

2/15 以降に解禁

chartr(intToUtf8(65:90),intToUtf8(c(66:90,65)),a) だったんだけど

24文字解 chartr("A-Z", "B-ZA", a) で十分だ(intToUtf8 を知ったばかりだったので,珍し物好きの気がむくむくと)

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

数学II・数学B 第3問 の出題ミス

2015年01月20日 | 雑感

自然数 n に対し,2n の 1 の位の数を an とする。また,数列 bn は b1 = 1,bn+1 = an bn / 4 (n=1,2,3,...) を満たすとする。

中略 このことから,すべての自然数 n に対して a = an となることがわかる。

オに当てはまるものを選べ

5n, 4n+1, n+3, n+4, n+5

正確な記載を参照のこと

正解は 5n でもよいし n+4 でもよい

> library(gmp)
> n = 1:500
> an = as.bigz(2)^n %% 10
> an[-(1:4)]
Big Integer ('bigz') object of length 499:
  [1] 2 4 8 6 2 4 8 6 2 4 8 6 2 4 8 6 2 4 8 6 2 4 8 6 2 4 8 6 2 4 8 6 2 4 8 6 2 4 8 6 2 4 8 6 2 4 8 6 2 4
     略
[478] 4 8 6 2 4 8 6 2 4 8 6 2 4 8 6 2 4 8 6
> an[seq(5,length(an), by=5)]
Big Integer ('bigz') object of length 100:
  [1] 2 4 8 6 2 4 8 6 2 4 8 6 2 4 8 6 2 4 8 6 2 4 8 6 2 4 8 6 2 4 8 6 2 4 8 6 2 4 8 6 2 4 8 6 2 4 8 6 2 4 8 6 2
 [54] 4 8 6 2 4 8 6 2 4 8 6 2 4 8 6 2 4 8 6 2 4 8 6 2 4 8 6 2 4 8 6 2 4 8 6 2 4 8 6 2 4 8 6 2 4 8 6

お粗末だなあ。受験生は最初の数個を実際に計算して,そこから規則性を導こうとする。そのほかの規則があるかどうかに気づく機会はない。

この短絡思考を逆用して不正解に導くという高等戦略もあるのだが,今回は,その戦略が失敗してしまった。

嘆かわしい...

でもねえ,こんな小細工(?)でふるい分けようという魂胆そのものがおかしいわけだ...

だからといって,総合的な学力に基づいて選抜しようというのはもっともっと難しい訳なんだが。

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

ヒストグラムの一部分に色をつける

2015年01月20日 | ブログラミング

いろいろ拡張した。

いつか必要になったときに備えて,メモしておく。

histcol2 = function(obj, lower, upper, probability = FALSE, ...) {
  if (class(obj) != "histogram") {
    obj = hist(data, plot = FALSE)
  }
  y = if (probability) obj$density else obj$counts
  breaks = obj$breaks
  n = length(y)
  paint = (upper[1] > breaks) & (breaks >= lower[1])
  if (length(lower) > 1) {
    for (i in 2:length(lower)) {
      paint = paint | ((upper[i] > breaks) & (breaks >= lower[i]))
    }
  }
  for (i in 1:n) {
    if (paint[i]) {
      rect(breaks[i], 0, breaks[i + 1], y[i], ...)
    }
  }
}
set.seed(123)
x = rnorm(1000)
a = hist(x)
histcol2(a, lower = -1, upper = 1, density = 20, col = "red")
a = hist(x, probability = TRUE, main = "")
histcol2(a, lower = -1, upper = 1, probability = TRUE, col = "aliceblue")
a = hist(x, breaks = c(-3, -2, seq(-1, 2, by = 0.5), 3, 4), probability = TRUE, main = "")
histcol2(a, lower = c(1, -3), upper = c(3, -2), probability = TRUE, col = "aliceblue")

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

「標本の大きさ」と「標本数」は違うよ!

2015年01月20日 | 統計学

> 今年のセンター試験の数学II・数学B第5問の「確率分布と統計的な推測」解説続き。(3)は設問文の日本語がややこしいのでしっかり読み解くことが重要。標本の数が4倍になっても、信頼区間の幅はその平方根倍(すなわちここでは2倍)にしかならないことに注意が必要。

「標本の数」という,専門用語にもない呼び方をしているが,「標本の大きさ(サンプルサイズ)」のことを言っているのは明らか。

こんな基本的なことが区別できない人が,統計学の問題の解説をしている...

追記:

あたりまえだが,問題文にはちゃんと「大きさ n の無作為標本を抽出する」って書いてある。

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

英語のリスニング機器不良

2015年01月20日 | 雑感

センター試験の英語(リスニング)において,全国で 86 人が機器の不良にぶち当たったのをメディアが「トラブル続出」と報道した件で,「わずか0.017%(86/520000*100)で,『続出か?』」とつぶやき,メディアを dis って悦に入っている。

こいつら馬鹿か?こういうものは,% で評価するもんじゃないよ。件数だ。統計学から言ってもあたりまえだ。

52 万人の客がハンバーガーを買って,その中で 86 人のハンバーガーにゴキブリが入っていたら,そのハンバーガー屋はつぶれるだろう。

食品業界などでは(どんな業界でも同じだろう),100 万個の製品中の不良品を 1 ~ 5 個に抑える(さらに,それら不良品は最終的には出荷されない)ような品質管理をしている。

 センター試験で使用される機器の不良率が 0.017% というのは,「高すぎる」のである。

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

再帰関数と非再帰関数

2015年01月19日 | ブログラミング

n 段の階段を,「1段,2段まとめ(1段飛ばし),3段まとめ(2段飛ばし)で昇る」を組み合わせると,昇り方は何通りあるか。

再帰的に考えればよいのだが,まじめにやると,n が大きくなると大変時間がかかる。

kaisuu = function(n) {
  if (n == 0) {
    return(1)
  } else if (n == 1) {
    return(Recall(n-1))
  } else if (n == 2) {
    return(Recall(n-1) + Recall(n-2))
  } else {
    return(Recall(n-1) + Recall(n-2) + Recall(n-3))
  }
}

> system.time(print(kaisuu(20)))
[1] 121415
   ユーザ   システム       経過  
     0.428      0.003      0.444
> system.time(print(kaisuu(30)))
[1] 53798080
   ユーザ   システム       経過  
   167.456      0.850    173.378

余計な再帰呼び出しが少なくなるようにしただけでも,実行時間は短くなる。

kaisuu2 = function(n) {
  if (n == 0) {
    return(1)
  } else if (n == 1) {
    return(1)
  } else if (n == 2) {
    return(2)
  } else {
    return(Recall(n-1) + Recall(n-2) + Recall(n-3))
  }
}

> system.time(print(kaisuu2(20)))
[1] 121415
   ユーザ   システム       経過  
     0.203      0.002      0.212
> system.time(print(kaisuu2(30)))
[1] 53798080
   ユーザ   システム       経過  
    90.573      0.597     95.541

しかし,再帰を用いないプログラムは問題なく速い。

kaisuu3 = function(n) {
  way = 0
  for (s1 in 0:n) {
    for (s2 in 0:(n%/%2)) {
      for (s3 in 0:(n%/%3)) {
        if (s1+2*s2+3*s3 == n) {
          way = way + factorial(s1+s2+s3)/factorial(s1)/factorial(s2)/factorial(s3)
        }
      }
    }
  }
  way
}

> system.time(print(kaisuu3(20)))
[1] 121415
   ユーザ   システム       経過  
     0.002      0.000      0.002
> system.time(print(kaisuu3(30)))
[1] 53798080
   ユーザ   システム       経過  
     0.014      0.000      0.023

しかし,まあ,この問題は「1段と2段で昇る」というのは Web ページにもよく見られるが,フィボナッチ数列の一般化ということで,以下のようにプログラムすれば所要時間は「検出限界以下」である。

kaisuu4 = function(n) {
  x = 1
  y = 2
  z = 4
  for (i in 4:n) {
    a = x+y+z
    x = y
    y = z
    z = a
  }
  z
}









コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

掛け算表

2015年01月16日 | ブログラミング

R で 9×9 掛け算表を表示する(ネタ)」で,

> こんなに短く書ける言語は他には存在しない!?

として,1:9%*%t(1:9) が書かれていた(空白を除いて12文字)が,

1:9%o%1:9 がもっと短い(同,9文字)

コメント (1)
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

1月,2月,...,12月

2015年01月15日 | ブログラミング

メモ

> format(ISOdate(2000, 1:12, 1), "%B")
 [1] "1月"  "2月"  "3月"  "4月"  "5月"  "6月"  "7月"  "8月"  "9月"  "10月" "11月" "12月"

> paste(1:12, "月", sep="")
 [1] "1月"  "2月"  "3月"  "4月"  "5月"  "6月"  "7月"  "8月"  "9月"  "10月" "11月" "12月"

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

規則性のある文字列出力(4)

2015年01月14日 | ブログラミング

規則性のある文字列出力(2)」だが,やっとのことで,37 文字解に到達できた。
二番煎じなので,投稿はしません。

締め切りが過ぎたので,解を提示する

cat(intToUtf8(65:90+32*(0:207%%8>0)))

rawToChar(as.raw(*)) が無駄だということで,intToUtf8(*) で十分ということだった

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

PVアクセスランキング にほんブログ村

PVアクセスランキング にほんブログ村