2016年12月8日木曜日

「エレガントな問題解決」演習問題 2.1.27(ひっかけ問題)(c)の解答

いつもは、↓別のブログの方に載せていたのですが、

「エレガントな問題解決」演習問題 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 × 4C3

ですね。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とか小数とか、そういうのを引数にするのは禁止ね)

じゃあ、いつ50%を超えるか、やっていきましょう。

> ProbMoreThan3(5)
[1] 0.5

> ProbMoreThan3(6)
[1] 0.65625

> ProbMoreThan3(7)
[1] 0.7734375

「確率が50%を超える」なので、正解は「6回」です。たぶん。

どのあたりが「ひっかけ」だったのでしょうか。

「3回以上」なのに、「ぴったし3回」の確率を求めしまうとか。

「50%を超える」なのに、50%になったとき(=5回)を答えてしまうとか。

それとも、何か別のひっかかる要素があって、私もまんまとそれにひっかかっているとか。
謎は深まります。

が、考えても分かる気がしないので、時間つぶしにグラフでも書いてみましょう。

n <- 3:20
p <- sapply(n, ProbMoreThan3)
plot(n,p)
abline(h=0.5, lty=2)


n回投げて3回以上表が出る確率

初めの方でぐっと上がって、それから1に漸近していくような感じですね。




2016年11月2日水曜日

Rで順列、組み合わせを計算する

最近知ったんですが、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という関数が使えます。1ずつ減っていく引数は、 ":" (コロン)を使えばよさそうですね。

ということで、

    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という関数があるので、そのまま使えます。

こちらも関数定義にしてみましょう。

# 組み合わせを計算する関数定義
Combi <- function(n, r){
  Permu(n, r)/factorial(r)
}


こちらも動作確認してみましょう。

> Combi(4, 0)
[1] 1
> Combi(4, 1)
[1] 4
> Combi(4, 2)
[1] 6
> Combi(4, 3)
[1] 4
> Combi(4, 4)

[1] 1

うん、よさそうです。




2016年6月16日木曜日

Rでデータフレームを部分抽出後、ベクトルや行列に変換する

例えば↓こんな感じで、各市町村の人口データがあったとします。

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

で、人口や男女比を比較するために、それぞれの市町村に対してグラフを描いてみようとしたとします。

じゃあ、1行ごとにループで回しながら、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"

構成要素としては、2つの数値の並びになっているもの、型がデータフレームなもんだから、barplotにそのまま渡すとエラーになってしまうのね。ちなみに、円グラフのpie関数なんかでも同じで、データフレームのままでは受け取ってくれない。

じゃあ、ベクトルに変換してあげよう。

barplot(as.vector(df[i, c("male", "female")], mode="numeric"))

as.vector後にbarplotに渡したもの

ラベル(maleとかfemale)の情報は消えちゃったけど、まあなんとかなったと。

エラーメッセージには「ベクトルか行列でなければなりません」とあったから、今度は行列にしてみよう↓


barplot( as.matrix(df[i, c("male", "female")]) )



as.matrix後にbarplotに渡したもの

列名が残っているおかげで、軸ラベルがそのまま出てくれるので、こっちの方がいいですね。
ちなみに↓こんな風に書いても、同じ動きになります。

barplot( t(t(df[1, c("male", "female")])) )


転置したときに行列に変換されるからですね。文字数的には少なく済みますが、あとで見たときや、他の人が見たときに混乱するかもしれないので、こういう書き方はやめた方がよさそうです。
で、ここからは補足になるんですが、もともとデータフレームを入力とするggplotを使ってみたらどうだろうかと。

前述のデータフレームはwideフォーマットなので、いったんlongフォーマットに変換します。

library(reshape2)

# longフォーマットに変換
df_l <- melt( df, id.vars="city", variable.name="gender",
              value.name="population" )
df_l
  city gender population
1  A市   male         95
2  B町   male         49
3  C村   male         26
4  D市   male         45
5  A市 female         99
6  B町 female         51
7  C村 female         24
8  D市 female         55

「dl_l」という名前のロングフォーマットのデータフレームに変換できたので、これをggplotに渡します。

library(ggplot2)


# A市の行だけ取り出してプロット

ggplot( df_l[city=="A市", ], aes(x=gender, y=population)) +
  geom_bar(stat="identity")



1行だけ取り出してggplot
うん、想定通りに表示できました。

じゃあ、これをループで回す? いやいや、ggplotなら、いっぺんにプロットできますよね。

ggplot( df_l, aes(x=city, y=population, fill=gender)) +
  geom_bar(stat="identity", position="dodge")



各市町村の男女人口や比率を比較したいなら、この可視化がよさげですね。

でも、このデフォルトのmaleとfemaleの色の割り当ては、無用な勘違いを生んでしまいそうです。
ggplotで、塗りつぶし(fill)の色を変更する方法については、

Rグラフィックスクックブック ―ggplot2によるグラフ作成のレシピ集
Rグラフィックスクックブック ―ggplot2によるグラフ作成のレシピ集

↑この書籍の「レシピ3.2 棒をグループ化する」のところなんかに、カラーパレットとしてbrewerパレットを使うレシピが載っていたりします。

その方法を使ってみると↓


はい、それっぽい色になりました。

やり方については、上記の書籍を買って読んでみてください。(最後はアフィリエイトかよ!)




2016年3月30日水曜日

Rを使って福岡市の人口密度ランキング(町丁目単位)

福岡県の市単位(政令指定都市は区単位)の人口密度単位のランキングは、↓こちらにあって、

福岡県の人口密度番付 - 都道府県・市区町村ランキング【日本・地域番付】

  1位 福岡市 中央区(11770人/㎢)
  2位 福岡市 城南区( 8031人/㎢)
  3位 福岡市 南区 ( 7976人/㎢)

となっております。トップは予想通りの福岡市中央区。

ただ、中央区と言ってもそれなりに広いし、人口密度にはムラがあるでしょうし、もしかしたら他の区にもギュギュっと人の密集しているエリアがあるかもしれません。

ということで、Rを使って町丁目単位(博多駅前1丁目とか、天神2丁目とか)でのランキングを出してみました。

とりあえず結果をどうぞ↓



中央区全体の人口密度は 11770人/㎢ でしたが、中央区荒戸2丁目で人口密度を計算すると 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 <-readShapePoly("merge.shp", proj4string=pj)
 
# 必要な項目を取り出しデータフレームにする
df <- data.frame(町丁名=shp$MOJI, 面積=shp$AREA, 人口=shp$JINKO)
 
# 人口密度の列を追加
df$人口密度 <- df$人口 / (df$面積 / 1000^2) # 1平方キロあたり
 
# ソートして、上位50だけ取り出す
df_sort <- df[order(df$人口密度, decreasing=T), ]
top50 <- df_sort[1:50, ]
 
# 画像ファイルに出力
png("人口密度ランキング.png", width=640, height=1280)
ggplot(top50, aes(x=reorder(町丁名, 人口密度), y=人口密度)) +
  geom_bar(stat="identity") +
  ggtitle("福岡市の人口密度TOP50") +
  ylab("人口密度 (人/㎢)") +
  theme( plot.title=element_text(size=20),
         axis.title.x=element_blank(),
         axis.title.y=element_blank(),
         axis.text.x=element_text(size=20),
         axis.text.y=element_text(size=20)  ) +
  geom_text(aes(label=paste(round(人口密度), "人/㎢"), y=人口密度-5000),
            size=5, colour="white") +
  coord_flip()
dev.off()