Rプログラミングの小ネタ

いつもは、↓別のブログの方に載せていたのですが、 「エレガントな問題解決」演習問題 2.1.25の解答（分母が3つの項の積になっている数列の和）: 主張 今回はRを使うので、こちらのブログで。 「エレガントな問題解決」 P29 2.1.27（ひっかけ問題）(c) 偏りのないコインを投げて、表が3回以上出る確率が50%を超えるには、最低何回投げればよいだろうか。 自分が題意を理解できているのか、いまいち自信がありません。私は、ひっかかっているのか？ とにかく解答してみます。 コインを1回投げて、表が3回以上出ることは、もちろんありません。 なので、考慮に値するのは3回以上投げるケースから。 3回投げて、表が3回出る確率をRで求めてみると、 > (1/2)^3 [1] 0.125 なので、12.5%。これは50%を超えていませんね。 では、4回投げる場合。表が3回出る確率と、表が4回出る確率を足せばいいですね。 4回投げて、表が3回出る確率は、（Cは組み合わせの記号だと思ってください） (1/2)^4 × 4 C 3 ですね。Rでは、組み合わせ数を求めるのにchooseコマンドが使えるので、 > (1/2)^4 * choose(4,3) [1] 0.25 4回投げて、表が3回出る確率は、 > (1/2)^4 * choose(4,4) [1] 0.0625 上記の2つを足しても、まだ50%は超えませんね。 いちいち足すのもまどろっこしいので、関数にしてみましょう。 # # n回投げたときに、3回以上でる確率を返す関数 # ProbMoreThan3 <- function(n){   s <- 0    # 組み合わせ数の合計値格納用   for(i in 3:n){     s <- s + choose(n,i)   }   p <- 1/2^n * s   return(p) } 試しに呼んでみましょう。 > ProbMoreThan3(3) [1] 0.125 > ProbMoreThan3(4) [1] 0.3125 よさそうです。（エラー処理は入れていないので、2とか小数とか、そういうのを引数にす

最近知ったんですが、Rで順列（Permutation）や、組み合わせ（Combination）を計算するデフォルトの関数ってないんですね。 【2016年11月18日追記 開始】 勘違いしていました。訂正します。 組み合せの数を計算するchooseという関数があります。 例えば、4個から2個選ぶ組み合わせの数は > choose(4,2) [1] 6 みたいな感じで求めることができます。 【2016年11月18日追記 終了】 そういうパッケージはあるみたいですけど、わざわざパッケージを使うまでもないような気もしないでもないです。 高校数学を復習する意味でも、自分で実装してみましょう。 まず、順列。 n個の中からk個を取り出して、順番を意識して並べていくような場合の数ですね。 nから初めて、1ずつ減らしながら、r個分かければ計算できます。     n × (n-1) × (n-2) × ... × {n-(r-2)} × {n-(r-1)} 上記をやりたい場合、受け取った引数を全部かけてくれるprodという関数が使えます。１ずつ減っていく引数は、 ":" （コロン）を使えばよさそうですね。 ということで、     prod(n:(n-r+1)) とやればいいのですが、これだと r が 0 のときにうまくいきません。0個を選び出す順列の数は 1 なので、関数にしてしまって、↓こんな感じで、どうでしょうか？ # 順列を計算する関数定義 Permu <- function(n, r){   ifelse( r==0, 1, prod(n:(n-r+1)) ) } 動作確認してみましょう。 > Permu(4, 0) [1] 1 > Permu(4, 1) [1] 4 > Permu(4, 2) [1] 12 > Permu(4, 3) [1] 24 > Permu(4, 4) [1] 24 うん、よさそうです。 次は、組み合わせ。 nからr個選ぶ順列を、rの階乗で割れば（r個の順番を無視する）、OKですね。 階乗はfactorialという関数があるので、そのまま使えます。 こちらも関数定義にしてみましょう。 # 組み

例えば↓こんな感じで、各市町村の人口データがあったとします。 city   <- c("A市", "B町", "C村", "D市") male   <- c(95, 49, 26, 45) female <- c(99, 51, 24, 55) df <- data.frame(city, male, female) df   city male female 1  A市   95     99 2  B町   49     51 3  C村   26     24 4  D市   45     55 で、人口や男女比を比較するために、それぞれの市町村に対してグラフを描いてみようとしたとします。 じゃあ、１行ごとにループで回しながら、barplotあたりを使って、グラフを描いてみようとしまして、↓こんな風にスクリプトを書くと（iがループのカウンタだと思ってね）、 # i行目の男女人口を棒グラフにしたい・・・ barplot(df[i, c("male", "female")])  barplot.default(df[i, c("male", "female")]) でエラー:    'height' はベクトルか行列でなければなりません なんでだ？と思って、切り出される部分をチェックしてみると、 df[i, c("male", "female")]   male female 1   95     99 class(df[i, c("male", "female")]) [1] "data.frame" 構成要素としては、２つの数値の並びになっているもの、型がデータフレームなもんだから、barplotにそのまま渡すとエラーになってしまうのね。ちなみに、円グラフのpie関数なんかでも同じで、データフレームのままでは受け取ってくれない。 じゃあ、ベクトルに変換してあげよう。 barplot( as.vector(

福岡県の市単位（政令指定都市は区単位）の人口密度単位のランキングは、↓こちらにあって、 福岡県の人口密度番付 - 都道府県・市区町村ランキング【日本・地域番付】 　　1位 福岡市 中央区（11770人/㎢） 　　2位 福岡市 城南区（ 8031人/㎢） 　　3位 福岡市 南区　（ 7976人/㎢） となっております。トップは予想通りの福岡市中央区。 ただ、中央区と言ってもそれなりに広いし、人口密度にはムラがあるでしょうし、もしかしたら他の区にもギュギュっと人の密集しているエリアがあるかもしれません。 ということで、Rを使って町丁目単位（博多駅前１丁目とか、天神２丁目とか）でのランキングを出してみました。 とりあえず結果をどうぞ↓ 中央区全体の人口密度は 11770人/㎢ でしたが、中央区荒戸２丁目で人口密度を計算すると 49294人/㎢ と、4倍以上の数字を叩き出しました。 また、上位5位の荒戸、平尾、薬院、高砂はいずれも中央区の町名ですが、6位に入っている愛宕浜は西区にある町です。 西区と言えば福岡市の中では辺境ですし（注：著者の個人的なイメージです）、前述の市区単位のランキングでは17位とかなり下位にランキングされてました、 実は、この愛宕浜2丁目は室見川をはさんで早良区のすぐ対岸にあって、比較的都市部に近い場所なんですよね。そして、ほとんどマンションしかありません（他には中学校とマルキョウがある）。大抵のマンションは10階以上あるので、おのずと人口密度が高くなるわけですね。 ↓Google Earthで見てみるとこんな感じ では、最後にRスクリプトを載せておきます。 シェープファイルは↓こちらに書いてある方法で入手して、 e-Stat(統計局)で公開されているShapeファイルを、Rで表示する - Rプログラミングの小ネタ ↓このやり方でマージしました。 Rで複数のシェープファイルを結合する - Rプログラミングの小ネタ library ( maptools ) library ( ggplot2 )   # シェープファイルを読み込む pj <- CRS ( "+proj=longlat +datum=WGS84" ) shp <-r