当選確率nのくじをn回ひいた時に少なくとも1回は当る確率の計算について
先日このツイートを見て、念のためそれを確かめてみた。
ガチャの予算計算をするときに考える覚えることは3個だけあればいいと考えてる。
— SCP-███ Lolida(ろりいだ) (@Lolida21) 2017年1月12日
ひとつはさっき言った
・n分の1のガチャをn回引いて当たりが出る確率は63%
で、それに追加して
・3n回引くと95%で当たりが出る
・5n回引くと99%で当たりが出る
だけ覚えれば破産はしないと思う
n回中、少なくとも1回は当る確率
「n回の試行中、少なくとも1回は当たりが出る確率」は、「n回の試行中、1回も当たらない確率」を「1（＝100％）」から引けばいい。 *1
当る確率が「p」だとして、1回の試行で当たらない確率は「1- 1/p」。 n回連続で当たらない確率は、「（1-1/p）n」（累乗計算）となります。
試しにPythonで数式計算をしてみる
Pythonで数式計算をするライブラリは「sympy」
from sympy import *
sympyでは先ず使用する変数を定義します。 ここでは1回の当選確率を「p」。求めたい値であるn回試行中1回は当る確率を「p1」。試行回数を「n」としました。
n,p,p1=symbols("n p p1")
先ほどの計算式を記述
p1 = 1 - (1-1/p)**n
また、今の設定では、当選確率の逆数の回数分だけくじを引くことになっているので「p == n」を代入し、それを新たに「p1_2」と定義します。
p1_2 = p1.subs(p,n)
print(p1_2)
すると、「1/n」の確率で当るくじを「n」回ひいた時に、少なくとも1回は当るが出る確率は
「-(1 - 1/n)**n + 1」
とわかりました。 先ほどの手計算の数式と同じです。
このp1_2に具体的に「n==2」なり「n==100」なりを代入した結果を出力すれば、求めたい数値が出てきます。 「evalf()」で分数が小数点表示になります
print(p1_2.subs(n,2).evalf()) print(p1_2.subs(n,6).evalf()) print(p1_2.subs(n,100).evalf())
当選確率1/2で2回だと「0.750000000000000」 1/6の確率で6回だと「0.665102023319616」 1/100の確率で100回だと「0.633967658726771」
値を求めるだけならRでできる
数式を求めることに意義が薄いと思えば最初からRを使えば手っ取り早いです。
当たり／外れ の2種類しかない試行は二項分布になります。 確率を求めるのは、「d～」で、二項分布は「binom」。 引数は、「当たりの確率」、「試行回数」、「当る確率」の3つで、先ほど同様、1回も当たらない（＝＝０）確率を1から引けば良いです。
1-dbinom(0,2,1/2) 1-dbinom(0,4,1/4) 1-dbinom(0,6,1/6) 1-dbinom(0,100,1/100)
> 1-dbinom(0,2,1/2) [1] 0.75 > 1-dbinom(0,4,1/4) [1] 0.6835938 > 1-dbinom(0,6,1/6) [1] 0.665102 > 1-dbinom(0,100,1/100) [1] 0.6339677
それはそうと、さっきと数字が違うじゃないですか？
さっきのツイートの数字は、63％と断言されていましたが、実際には違う、、、でも1/100で100回の時はとても近い
Rで繰り返し処理
というわけで、n==2 から n==100 までの夫々の確率をデータフレームに格納して、プロットしてみます。
繰り返し処理でデータフレームに追記していく場合、rbindを使います。
df <- data.frame() for (n in 2:100){ dat <- data.frame(1- dbinom(0,n,1/n)) df <- rbind(df,dat) } colnames(df) <- c("prob") df$x <- rep(2:100)
結果を出力すると、試行回数が上がるに連れて一定の数字に修練しているようです。
> head(df) prob x 1 0.7500000 2 2 0.7037037 3 3 0.6835938 4 4 0.6723200 5 5 0.6651020 6 6 0.6600833 7 > tail(df) prob x 94 0.6340653 95 95 0.6340450 96 96 0.6340250 97 97 0.6340055 98 98 0.6339864 99 99 0.6339677 100
グラフにしてみる
ggpotでグラフ化
p <- ggplot(data=df,aes(x=x,y=prob)) p <- p + geom_bar(stat="identity") p <- p + ggtitle("当たりの確率の逆数と同じ回数だけひいたら１回は辺りが出る確率") p <- p + geom_text(aes(x=100, y=df[df$x==100,]$prob, label=round(df[df$x==100,]$prob*100,2)), vjust=-1) plot(p) ggsave(plot=p,file="20170115.png",h=4,w=6,scale=2)
ほぼほぼ63.4%に収斂しています。
それでは、「3n回」や「5n回」は？
繰り返し計算中の確率計算部分で試行回数を「n3」や「n5」にするだけです
dat <- data.frame(1- dbinom(0,n*3,1/n))
3倍の回数なら
5倍の回数なら
実際のところ、確率1％であっても100回も引く人はいないと思いますが、、、
これだけ回数を増やせばポアソン分布である
ガチャの確率が一般的にどれくらいなのか知りませんが、 二項分布は確率が低く、試行回数が増えればポアソン分布で近似できる、ということがわかっています。
例えば、当選確率pが「1/100」で試行回数nが「100回」なのであれば、
lambda = n * p は、「100*(1/100)」
同様に「300回」なら「100(1/300)」、「500回」なら「500(1/100)」を引数として
1 - dpois(0,(100*(1/100))) 1 - dpois(0,(300*(1/100))) 1 - dpois(0,(500*(1/100)))
とすれば、先ほど収束していった値が出てきます。
> 1 - dpois(0,(100*(1/100))) [1] 0.6321206 > 1 - dpois(0,(300*(1/100))) [1] 0.9502129 > 1 - dpois(0,(500*(1/100))) [1] 0.9932621
コード全体
from sympy import * n,p,p1=symbols("n p p1") p1 = 1 - (1-1/p)**n p1_2 = p1.subs(p,n) print(p1_2) print(p1_ans.evalf()) print(p1_2.subs(n,6).evalf()) print(p1_2.subs(n,100).evalf())
library("dplyr", lib.loc="C:/hogehoge/r/R-3.3.2/library") library("reshape2", lib.loc="C:/hogehoge/r/R-3.3.2/library") library("ggplot2", lib.loc="C:/hogehoge/r/R-3.3.2/library") 1-dbinom(0,2,1/2) 1-dbinom(0,4,1/4) 1-dbinom(0,6,1/6) 1-dbinom(0,100,1/100) df <- data.frame() for (n in 2:100){ dat <- data.frame(1- dbinom(0,n,1/n)) df <- rbind(df,dat) } colnames(df) <- c("prob") df$x <- rep(2:100) head(df) tail(df) p <- ggplot(data=df,aes(x=x,y=prob)) p <- p + geom_bar(stat="identity") p <- p + ggtitle("当たりの確率の逆数と同じ回数だけひいたら１回は辺りが出る確率") p <- p + geom_text(aes(x=100, y=df[df$x==100,]$prob, label=round(df[df$x==100,]$prob*100,2)), vjust=-1) plot(p) ggsave(plot=p,file="20170115.png",h=4,w=6,scale=1) df_3 <- data.frame() for (n in 2:100){ dat <- data.frame(1- dbinom(0,n*3,1/n)) df_3 <- rbind(df_3,dat) } colnames(df_3) <- c("prob") df_3$x <- rep(2:100) print(df_3) p3 <- ggplot(data=df_3,aes(x=x,y=prob)) p3 <- p3 + geom_bar(stat="identity") p3 <- p3 + ggtitle("当たりの確率の逆数との３倍の回数をひいたら１回は辺りが出る確率") p3 <- p3 + geom_text(aes(x=100, y=df_3[df_3$x==100,]$prob, label=round(df_3[df_3$x==100,]$prob*100,2)), vjust=-1) plot(p3) ggsave(plot=p3,file="20170115_3.png",h=4,w=6,scale=1) df_5 <- data.frame() for (n in 2:100){ dat <- data.frame(1- dbinom(0,n*5,1/n)) df_5 <- rbind(df_5,dat) } colnames(df_5) <- c("prob") df_5$x <- rep(2:100) print(df_5) p5 <- ggplot(data=df_5,aes(x=x,y=prob)) p5 <- p5 + geom_bar(stat="identity") p5 <- p5 + ggtitle("当たりの確率の逆数との３倍の回数をひいたら１回は辺りが出る確率") p5 <- p5 + geom_text(aes(x=100, y=df_5[df_5$x==100,]$prob, label=round(df_5[df_5$x==100,]$prob*100,2)), vjust=-0.5) plot(p5) ggsave(plot=p5,file="20170115_5.png",h=4,w=6,scale=1) 1 - dpois(0,(100*(1/100))) 1 - dpois(0,(300*(1/100))) 1 - dpois(0,(500*(1/100)))
*1:数学が面白いのはこの手の逆転の発想が普通に出てくることですよね
日本の労働者の低い生産性_再現可能なものへの転換
年末年始最初に読んだ本は、
![デービッド・アトキンソン 新・所得倍増論 潜在能力を活かせない「日本病」の正体と処方箋 [ デービッド・アトキンソン ]](https://thumbnail.image.rakuten.co.jp/@0_mall/book/cabinet/6353/9784492396353.jpg?_ex=128x128 "デービッド・アトキンソン 新・所得倍増論 潜在能力を活かせない「日本病」の正体と処方箋 [ デービッド・アトキンソン ]")
デービッド・アトキンソン 新・所得倍増論 潜在能力を活かせない「日本病」の正体と処方箋 [ デービッド・アトキンソン ]
- ジャンル: 本・雑誌・コミック > ビジネス・経済・就職 > 経済・財政 > 日本経済
- ショップ: 楽天ブックス
- 価格: 1,650円
でした。
昨今では日本の高度経済成長の要員はモーレツ社員の頑張りだけでなく、「円の固定相場制」「人口ボーナス」などの影響が大きかった、というのは通説になりつつあるのではないかと思いますが、本書は「人口ボーナス」の要員の大きさを自覚すべきだという主張でした。
そしてこのまま生産性が低いママでは各国に置いてけぼりになります。
では、どうしたら生産性は上がるのか？
デービッド アトキンソン氏の本は、アナリストらしく、具体的な改善方法というよりも、マネジメント層が改革をせざるを得なくなるように株式市場のプレッシャーを与えるように構造を変えよう、という主張で、その先は経営者が自分で考えることになっています。
普段、日経CPの元編集長 木村さんのツイートをよく見ていますが、日本は”失われた20年”の間に随分と効率化から置いてけぼりになっているのだと毎朝思い知らされる事例をツイートしてくれます。
ITで業務改革と言うけれど、本当に必要なのはマネジメント改革。ITを活用すれば、単なるオーバーヘッドの管理職を一掃できるし、効果は抜群のはず。ただ、抵抗勢力管理職だし、その親玉は社長。だから、絶望的に難しい。唯一可能なのは、会社が潰れかけて外資に買収された時ぐらいだな。
— 木村岳史(東葛人) (@toukatsujin) 2017年1月3日
カスタマイズするのは、悪
何年か前に読んでなるほど！と思った本に
![V字回復の経営 2年で会社を変えられますか （日経ビジネス人文庫） [ 三枝匡 ]](https://thumbnail.image.rakuten.co.jp/@0_mall/book/cabinet/5321/53219342.jpg?_ex=128x128 "V字回復の経営 2年で会社を変えられますか （日経ビジネス人文庫） [ 三枝匡 ]")
があります。
小説ですがとてもためになります。
本書では効率化を図るために、部品のカスタマイズ＝「すり合わせ」を止めて、共通して売れるものを生産し、それを買ってくれる客を探すのだ、という主張をしていました。
「すり合わせ」とは逆方向に進むインダストリー4.0
中間部品の標準化の重要性は世界的に肯定されている流れです。
![日本版 インダストリー4.X 日本の強みを活かすIoT革命 [ 奥 雅春 ]](https://thumbnail.image.rakuten.co.jp/@0_mall/book/cabinet/4489/9784478084489.jpg?_ex=128x128 "日本版 インダストリー4.X 日本の強みを活かすIoT革命 [ 奥 雅春 ]")
にも、ドイツのインダストリー4.0は工場の自動化などではなく、中小企業が納品する中間生産財の規格を標準化していることの方が重要であるという指摘がありました。
再利用可能という考え方
規格化や標準化による効率化は、再利用可能と言い換えても良いと思います。
一つ一つ、個別対応してすり合わせをした部品は他の顧客には使えませんが、標準化した部品は他の顧客にも使えます。
プログラムの世界
例えば、統計のレポーティングをExcelではなくRなどプログラムで行うのは、コードが再利用可能だからです。
お酒の世界
職人芸の世界でも獺祭のような事例がありました。職人の技を再利用可能にした取組です。
広告デザインにおける再利用可能性について
では広告の世界ではどのような分野で再利用可能なものへの変革が可能でしょうか？
例えば、デザインの世界では、一文字一文字字間を詰めるだとか、手作業で調整しているところが大変沢山あります。
一方でAdobeは最先端の技術を道入してどんどん人間の手間を省いています。
因みについ最近、人間の切抜き作業をしたのですが一昔前とは全然切抜きの精度が向上していて、小さいサイズであれば全くのおまかせでも大丈夫そうなくらいでした。
手間がかかっている方が良いという信仰
この人のツイートもよく読んでいますが、例えばこういうことなのだと思います。
一つ一つ特別仕様にすることが良いことだと思っている人は多いのではないかと思います。
「これは手間がかかっている！」ということを良しとするのをやめるべきでしょう。
手間がかかっていないことは消費者のメリットになる
生産ロットが多いということは、イニシャルコストが薄く分散されているため、同じ値段のスマートフォンでもロット数が多いものほど、良い機能が突っ込まれているわけです。
China eyes breakthroughs in SOE reform - Business - Chinadaily.com.cn
工芸品と工業製品は異なるものですので、手間がかかっていることが良い場合と悪い場合とがあるということだと思います。
参考
固定相場制の影響の重要性について
![経済は世界史から学べ！ [ 茂木誠 ]](https://thumbnail.image.rakuten.co.jp/@0_mall/book/cabinet/3648/9784478023648.jpg?_ex=128x128 "経済は世界史から学べ！ [ 茂木誠 ]")
- ジャンル: 本・雑誌・コミック > ビジネス・経済・就職 > その他
- ショップ: 楽天ブックス
- 価格: 1,650円
人口ボーナスについて
![移行期的混乱 経済成長神話の終わり （ちくま文庫） [ 平川克美 ]](https://thumbnail.image.rakuten.co.jp/@0_mall/book/cabinet/0250/9784480430250.jpg?_ex=128x128 "移行期的混乱 経済成長神話の終わり （ちくま文庫） [ 平川克美 ]")
デザインの職人芸について
![なるほどデザイン 目で見て楽しむデザインの本。 [ 筒井美希 ]](https://thumbnail.image.rakuten.co.jp/@0_mall/book/cabinet/5174/9784844365174.jpg?_ex=128x128 "なるほどデザイン 目で見て楽しむデザインの本。 [ 筒井美希 ]")
なるほどデザイン 目で見て楽しむデザインの本。 [ 筒井美希 ]
- ジャンル: 本・雑誌・コミック > ホビー・スポーツ・美術 > 美術 > デザイン
- ショップ: 楽天ブックス
- 価格: 2,200円
【備忘録】Rで矢印のグラフを描く
グラフにこんな矢印を入れたい時ってありますよね？
これは、因子分析をした時のプロットをイメージしています。
色々と調べた結果をデータの準備からプロットまで手順を追ってメモをしておきます。
データの準備
それっぽいデータを乱数で用意
取り敢えず正規分布の乱数を作成。
平均を0、分散を0.5として7つずつ
set.seed(123) mr1 <- rnorm(7,mean = 0,sd=0.5) mr2 <- rnorm(7,mean = 0,sd=0.5)
それっぽいイメージ項目の文字列を作成
names <- c("美しい","賢い","派手な","高貴な","個性的な","明るい","温かい")
データフレームを用意してdplyrのbind_colsを使ってデータフレームを作ります （別にcbindでも同じことです）
res <- NULL res <- as.data.frame(res) library("dplyr", lib.loc="C:/hogehoge/R-3.2.3/library") res <- bind_cols(as.data.frame(mr1),as.data.frame(mr2)) res$names <- names
結果、このようなデータを用意しました
mr1 mr2 names 1 -0.28023782 -0.63253062 美しい 2 -0.11508874 -0.34342643 賢い 3 0.77935416 -0.22283099 派手な 4 0.03525420 0.61204090 高貴な 5 0.06464387 0.17990691 個性的な 6 0.85753249 0.20038573 明るい 7 0.23045810 0.05534136 温かい
矢印をひく為のデータを用意
このサイトを参考にしました。
5 functions to do Principal Components Analysis in R · Gaston Sanchez
矢印には出発点とゴール地点とがあります。その値を持ったデータフレームを別途用意すれば良いのです。
arrowsという名称でデータフレームを用意します。
x1,x2が出発点なので"0"。今回のデータのレコード数をnrowで確認して、必要な数だけの0を並べます。 y1,y2はゴール地点なので描画するデータの値と同じです
arrows <- NULL arrows <- as.data.frame(arrows) x1 <- rep(0,nrow(res)) y1 <- rep(0,nrow(res)) arrows <- bind_cols(as.data.frame(x1),as.data.frame(y1)) x2 <- res$mr1 y2 <- res$mr2 arrows$x2 <- x2 arrows$y2 <- y2
こんなデータになります
x1 y1 x2 y2 1 0 0 -0.28023782 -0.63253062 2 0 0 -0.11508874 -0.34342643 3 0 0 0.77935416 -0.22283099
ggplotで描く
ggplotの描画範囲を大きくする
グラフの端で文字が消えないよう、x軸 y軸の最大値・最小値を少し大きめに取ります。
有効数字3桁として、0.1だけ絶対値を大きくとりました。
a <- round(max(res$mr1,res$mr2),3) + 0.1 b <- round(min(res$mr1,res$mr2),3) - 0.1
ベースの値をプロット
先ずはデータの値（今回は因子分析の因子負荷量をイメージしています。）をgeom_pointで、質問項目名をgeom_textで、先ほど決めたグラフの軸の値でxlim,ylimで指定します
g <- ggplot(data=res,aes(x=res$mr1,y=res$mr2)) g <- g + geom_point() g <- g + geom_text(aes(x=res$mr1,y=res$mr2),label=res$names,vjust=-1) g <- g + xlim(b,a) g <- g + ylim(b,a) plot(g)
するとこんなグラフになりますが、ちょっと見辛いです。
縦横の比率を同じにする
グラフの描画エリアが正方形であって欲しい場合があります。今回は、縦も横も因子負荷量ですから、同じスケールでないと正確にイメージできない可視化となってしまいます。
こちらのサイトを参考にしました。
ggplot2: きれいなグラフを簡単に合理的に - Watal M. Iwasaki
g <- g + coord_fixed(ratio=1)
を加えます。
また、負の値を取りますので、0の場所にx軸、y軸の軸がありません。 そこで、0の水平線と垂線を引いて見やすくします。色もブルーにしてみました。
g <- g + geom_hline(yintercept = 0,color="blue") g <- g + geom_vline(xintercept = 0,color="blue")
これでかなり見やすくなりました。
いよいよ矢印を引きます
こちらのサイトを参考にしました。
ggplot2 Quick Reference: geom_segment | Software and Programmer Efficiency Research Group
geom_segmentという関数を使って、先ほど作った矢印用のデータフレームarrowsを対象データとして使います。arrow=arrow()で直線が矢印に変わります。
今回は太さを1にし、色をグレーにしてみました。
g <- g + geom_segment(data=arrows,aes(x=x1,y=y1,xend=x2,yend=y2),colour="gray65",arrow=arrow(),size=1)
細かいこだわり
よくみると矢印が青い線の上に乗っかって青線が消えています。ggplotはグラフのデータにどんどんレイヤーを上乗せしていくので、上のレイヤーにしたデータは後で加えないと行けません。 ですので、geom_vlineとgeom_hlineは後ろに順番を入れ替えるべきです。
また、x軸、y軸の目盛りも黒文字の方が読み易いでしょう。
その結果のグラフが最初のグラフのです
全てのコード
set.seed(123) mr1 <- rnorm(7,mean = 0,sd=0.5) mr2 <- rnorm(7,mean = 0,sd=0.5) res <- NULL res <- as.data.frame(res) library("dplyr", lib.loc="C:/hogehoge/R-3.2.3/library") res <- bind_cols(as.data.frame(mr1),as.data.frame(mr2)) res$names <- c("美しい","賢い","派手な","高貴な","個性的な","明るい","温かい") arrows <- NULL arrows <- as.data.frame(arrows) x1 <- rep(0,nrow(res)) y1 <- rep(0,nrow(res)) arrows <- bind_cols(as.data.frame(x1),as.data.frame(y1)) x2 <- res$mr1 y2 <- res$mr2 arrows$x2 <- x2 arrows$y2 <- y2 library("ggplot2", lib.loc="C:/hogehoge/R-3.2.3/library") a <- round(max(res$mr1,res$mr2),3) + 0.1 b <- round(min(res$mr1,res$mr2),3) - 0.1 g <- ggplot(data=res,aes(x=res$mr1,y=res$mr2)) g <- g + geom_point() g <- g + geom_text(aes(x=res$mr1,y=res$mr2),label=res$names,vjust=-1) g <- g + xlim(b,a) g <- g + ylim(b,a) g <- g + coord_fixed(ratio=1) g <- g + geom_segment(data=arrows,aes(x=x1,y=y1,xend=x2,yend=y2),colour="gray65",arrow=arrow(),size=1) g <- g + geom_hline(yintercept = 0,color="blue") g <- g + geom_vline(xintercept = 0,color="blue") g <- g + theme(axis.text.x=element_text(angle=0,colour="black",size=12,hjust=1)) g <- g + theme(axis.text.y=element_text(angle=0,colour="black",size=12)) plot(g) ggsave(plot=g,file="20160830.png",dpi=300,width=4,height=1.5,scale=3)
【備忘録】交互作用プロットをggplot2で描く
今読んでいる本に「交互作用プロット」というものが出てきました。
マーケティング・データ分析の基礎 (シリーズ Useful R 3)
標準で入っている関数を使っていたので、せっかくならggplot2で描いてみようと思いました。
今回は、そのメモになります。
交互作用プロットとは？
交互作用プロットの例
2つの因子が存在した時に、片方の因子をx軸で、もう片方の因子をグラフの種類（色分け、点線／棒線等）で分けるというものです。
y軸は、平均値をプロットするそうです。
今回始めてみましたが、慣れないとちょっと見辛いですね。
交互作用について
見方や説明の詳細はこちらを参考にして下さい。
- 折れ線が並行かどうか ⇒ グラフの種類で分けた因子に効果がないか？あるか？
- 左右で上下があるかどうか ⇒ ｘ軸で分けた因子に効果があるかないか？
ということになります。
標準のグラフで描く
データの準備
このグラフはRのデータセットの中の「ギニアピッグの歯の成長に対するビタミン C の効果」のデータを使っています
len supp dose 1 4.2 VC 0.5 2 11.5 VC 0.5 3 7.3 VC 0.5
「ビタミン C の投与量」(dose)は、(0.5mg,1mg, 2mg)の3種類ありますので、今回の作図の都合により1mgだけ除きました
library("dplyr", lib.loc="C:/hogehoge/soft/r/R-3.2.3/library") df <- ToothGrowth %>% filter(dose==c(0.5,2))
また、「ビタミン C の投与量」(dose)は数値データなので因子（factor）に変換しておきます
df$dose <- as.factor(df$dose) str(df) 'data.frame': 20 obs. of 3 variables: $ len : num 4.2 7.3 6.4 11.2 5.2 18.5 25.5 32.5 21.5 29.5 ... $ supp: Factor w/ 2 levels "OJ","VC": 2 2 2 2 2 2 2 2 2 2 ... $ dose: Factor w/ 2 levels "0.5","2": 1 1 1 1 1 2 2 2 2 2 ...
標準の関数で描く
先ほどの本で紹介されていたinteraction.plot関数を使うと、先ほどのグラフになります。
interaction.plot(df$dose,df$supp,df$len,xlab="ビタミン C の投与量",trace.label="摂取法",ylab="造歯細胞(歯)の成長量")
ggplot2でやってみる
せっかくなら美しいggplot2で描いてみたいものですね。
データの準備
ggplot2の場合は、まずデータの集計が必要です。平均値をプロットするのですから、dplyr/reshape2で集計は簡単です。
先ずはデータを一度分解し、
library("reshape2", lib.loc="C:/hogehoge/r/R-3.2.3/library") df_m <- melt(df,id.vars = c("supp","dose")) head(df_m) supp dose variable value 1 VC 0.5 len 4.2 2 VC 0.5 len 7.3 3 VC 0.5 len 6.4
平均値に集計しなおします。dcastの際にmeanを使います。
df2 <- dcast(df_m,supp+dose~variable,mean) (df2) supp dose len 1 OJ 0.5 14.40 2 OJ 2 26.42 3 VC 0.5 6.86 4 VC 2 25.50
確かにこの表だけではちょっとイメージがつきにくいのかもしれません
ggplot2で描く
なお、ggplot2は非常に優秀なので、連続データでもないものを棒グラフで結ぶことを良しとしません。 そのためそのままgeom_lineと書くだけではエラーになってしまいます。
geom_path: Each group consists of only one observation. Do you need to adjust the group
aesthetic?
なので、無理矢理「supp」（摂取法）をグループだと認識させます。
g <- ggplot(df2,aes(x=dose,y=len,color=supp,group=supp)) g <- g + geom_line(size=2) g <- g + theme(axis.text.x=element_text(angle=0,colour="black",size=12,hjust=1)) g <- g + theme(axis.text.y=element_text(angle=0,colour="black",size=16)) plot(g) ggsave(plot=g,file="20160821_2.png",dpi=300,width=4,height=3,scale=2)
するとこのようなグラフになりました。
標準の関数はもっと長体がかかっているようなので、印象を併せる為にもう少し横長にしてみました。
ggsave(plot=g,file="20160821_3.png",dpi=300,width=4,height=1.5,scale=2)
分散分析もしておく
先ほどのデータで交互作用の有無を分散分析で確認しておきます
res <- aov(len~supp+dose+supp:dose,data=df) summary(res) Df Sum Sq Mean Sq F value Pr(>F) supp 1 89.5 89.5 5.748 0.0291 * dose 1 1175.0 1175.0 75.493 1.87e-07 *** supp:dose 1 54.8 54.8 3.519 0.0790 . Residuals 16 249.0 15.6 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
supp=摂取法 及び dose=ビタミンＣの投下量 については、5％の有意水準で”同じである＝効果がない”という帰無仮説は棄却されます。
しかし、supp:dose=摂取法と投下量の交互作用 については、5%以上なので、”同じである＝効果がない”ということは、”滅多にないことではない”ので、 交互作用があるとは言い難いようです。
データセットについて
余談ですが、こういう練習用のデータについてなかなか手元にドンピシャなものがないものです。
こういうまとめは本当に助かります
箱ひげ図とバイオリンプロットの印象の違い
データの可視化に於いて重宝されるグラフの一つに「箱ひげ図」というものがあります。
データの四分位点、中央値、最大値、最小値を記載するのでデータの分布が一覧できるというものです。
しかし箱ひげ図には重大な問題があります。
先ずはプロットしてみる
Rがあれば誰でも持っているアヤメのデータを使います。 ”がく”の長さを取り出して箱ひげ図にしてみます。
当然、最大値と四分位点、最小値と四分位点の間にもデータはあるのですが、箱が四分位点で終わってしまっているため、あたかもそこにはデータが無いような印象ではないでしょうか？
そこでバイオリンプロット
そんな箱ひげ図よりも便利な図としてバイオリンプロットというものが存在します。
ヒストグラムを縦にしたようなもので、これであればデータの凡その分布もわかります
箱ひげ図と比較してみる
こちらのサイトを参考にさせて頂きました。
バイオリンプロットの中にボックスプロットを入れると、より情報量が増えますね。
今回のコード
データの集計まで
ggplotで描けるようにreshapeパッケージでデータを変更
library("reshape2", lib.loc="C:/hogehoge/r/R-3.2.3/library") iris_m <- melt(iris,id.vars = c("Species"))
こんな感じになります
Species variable value 1 setosa Sepal.Length 5.1 2 setosa Sepal.Length 4.9 3 setosa Sepal.Length 4.7 4 setosa Sepal.Length 4.6 5 setosa Sepal.Length 5.0 6 setosa Sepal.Length 5.4
そのうち”がく”（＝Sepal）のデータだけをdplyrのfilter関数で抜き出します
library("dplyr", lib.loc="C:/hogehoge/r/R-3.2.3/library") iris_m_sl <- iris_m %>% filter(variable==c("Sepal.Length"))
tailで後半のデータを確認するとSepal.Lengthだけになっているのがわかります
Species variable value 145 virginica Sepal.Length 6.7 146 virginica Sepal.Length 6.7 147 virginica Sepal.Length 6.3 148 virginica Sepal.Length 6.5 149 virginica Sepal.Length 6.2 150 virginica Sepal.Length 5.9
ボックスプロット
ggplot2のgeom_boxplotを使います。 また先ほどのサイトを参考にして、平均値もサークルでプロットしておきます。
library("ggplot2", lib.loc="C:/hogehoge/r/R-3.2.3/library") iris_m_sl_g2 <- ggplot(iris_m_sl,aes(x=Species,y=value,color=Species,fill=Species)) iris_m_sl_g2 <- iris_m_sl_g2 + geom_boxplot(color="black") iris_m_sl_g2 <- iris_m_sl_g2 + stat_summary(fun.y=mean,geom = "point", fill="white",shape=21,size=3) iris_m_sl_g2 <- iris_m_sl_g2 + theme(axis.text.x=element_text(angle=0,colour="black",size=12,hjust=1)) iris_m_sl_g2 <- iris_m_sl_g2 + theme(axis.text.y=element_text(angle=0,colour="black",size=16)) plot(iris_m_sl_g2) ggsave(plot=iris_m_sl_g2,file="iris_m_sl_g_1.png",dpi=300,width=4,height=3,scale=2)
バイオリンプロット
バイオリンプロットはgeom_violinです。
デフォルトでは横幅が調整されてしまうので、データ量も比較したい場合は、scale="count"としておくと良いようです。
library("ggplot2", lib.loc="C:/hogehoge/r/R-3.2.3/library") iris_m_sl_g <- ggplot(iris_m_sl,aes(x=Species,y=value,color=Species,fill=Species)) iris_m_sl_g <- iris_m_sl_g + geom_violin(scale="count") iris_m_sl_g <- iris_m_sl_g + theme(axis.text.x=element_text(angle=0,colour="black",size=12,hjust=1)) iris_m_sl_g <- iris_m_sl_g + theme(axis.text.y=element_text(angle=0,colour="black",size=16)) plot(iris_m_sl_g) ggsave(plot=iris_m_sl_g,file="iris_m_sl_g_2.png",dpi=300,width=4,height=3,scale=2)
みんなのR
バイオリンプロットの方が情報量が多いよ！という話は、
![みんなのR 第2版 [ Jared P. Lander ]](https://thumbnail.image.rakuten.co.jp/@0_mall/book/cabinet/2159/9784839962159.jpg?_ex=128x128 "みんなのR 第2版 [ Jared P. Lander ]")
で指摘されていました。
箱ひげ図が先に広がってしまったためにバイオリンプロットは余りメジャーではないようですが（見た目も少し気持ち悪いですし） 誤った印象を与えないという意味ではもっと使って欲しいと思います。
RでGoogleAdWordsのCSVを読み込む時_文字コード_区切り文字_を指定してファイル読込
近年は誰でも安価にデータ分析ができるツールを手にする時代となりましたが、実務に使ってみようと思っても教科書通りにいかない、、ということはよくある現象かと思います。
リテラシーが高い人にとっては呼吸をするように無意識的に対応できるところが、私のような初心者にとっては躓きとなります。
今回は、AdWordsのレポートをCSVファイルで取得した時に文字化けで困ったので、その対応方法についてまとめてみました。
先ずは普通に読み込んでみた
コード
hogehoge <- read.csv("hogehoge.csv") head(hogehoge)
結果
Error in make.names(col.names, unique = TRUE) : invalid multibyte string at '<ff><fe><83>^JT<ec>0ﾝ0<fc>0ﾈ0' In addition: Warning messages: 1: In read.table(file = file, header = header, sep = sep, quote = quote, : line 1 appears to contain embedded nulls
エラーです
文字コードを確認する
Rではデフォルトの文字コード以外の文字コードのファイルを読み込む際には文字コードを指定して読み込む必要があります。
文字コードの確認には、nkfが便利です。 使い方やインストールはここでは詳解しませんので各自ご準備下さい。
コマンドプロンプトでファイルのある箇所まで行き、nkfで確認
結果
文字コードを指定して読込
read.tableでもread.csvでも文字コードを指定して読込ができます。「fileEncoding = "UTF-16"」を追加します。
コード
hogehoge <- read.csv("hogehoge.csv",fileEncoding = "UTF-16") head(hogehoge)
結果
広告レポート..2016.06.15.2016.07.14. 1 日\t広告の状態\t広告\t広告文 1\t広告文 2\t表示URL\tリンク先URL\t最終ページ URL\tモバイルの最終ページ URL\tトラッキング テンプレート\tカスタム パラメータ\tHeadline 1\tHeadline 2\t説明\tPath 1\tPath 2\tキャンペーン\t広告グループ\tステータス\tラベル\tキャンペーン タイプ\tキャンペーンのサブタイプ\tクリック数\t表示回数\tクリック率\t平均クリック単価\t費用\t平均掲載順位
文字化は解消されましたが、なんか変なデータになってしまいました。
先ずは、どうやら一行目にデータ以外の余分な行が入っているようです。人は何故データ以外の情報をCSVファイルに書いてしまうのでしょうか！！
冒頭のn行をスキップして読込
冒頭の何行かに余分な情報を入れていることはよくある話なので、スキップして読み込むという便利なオプション「skip=(行数)」があります。
コード
hogehoge <- read.csv("hogehoge.csv",fileEncoding = "UTF-16",skip=1) head(hogehoge)
結果*1
日.広告の状態.広告.広告文.1.広告文.2.表示URL.リンク先URL.最終ページ.URL.モバイルの最終ページ.URL.トラッキング.テンプレート.カスタム.パラメータ.Headline.1.Headline.2.説明.Path.1.Path.2.キャンペーン.広告グループ.ステータス.ラベル.キャンペーン.タイプ.キャンペーンのサブタイプ.クリック数.表示回数.クリック率.平均クリック単価.費用.平均掲載順位 ... 1 2016-07-13\tアクティブ\t広告名: 320x100_hogehoge.png; 320 x 100\t\t\t\t\t\t\t\t\t\t\t\t承認済み\t --\tディスプレイ ネットワークのみ\tすべての機能
まだなんだかよくわからない長い行として認識されています。
区切り文字を指定して読込
よく見ると「\t」という文字があります。正規表現でよく見る「タブ」のことだと推測されます。 CSVってカンマ区切りって意味じゃなかったっけ？と思いながら 区切り文字をタブに指定して 読込ます。
コード
hogehoge <- read.csv("hogehoge.csv",fileEncoding = "UTF-16",skip=1,sep="\t") head(hogehoge)
結果の前の方*2
日 広告の状態 広告 広告文.1 広告文.2 表示URL 1 2016-07-13 アクティブ 広告名: 320x100_hogehoge.png; 320 x 100 URL 2 2016-07-06 アクティブ 広告名: 468x60_hogehogec.png; 468 x 60 URL
後ろのほう*3
キャンペーン.タイプ キャンペーンのサブタイプ クリック数 表示回数 クリック率 1 ディスプレイ ネットワークのみ すべての機能 10 10000 __% 2 ディスプレイ ネットワークのみ すべての機能 300 100000 __%
どうやら１行目がヘッダー情報で、２行目からがデータのようです。
ヘッダーを指定する
既にread.csvが勝手にヘッダーを認識してくれているようですが、通常はヘッダーがあるよと指定して読込ます。
コード
hogehoge <- read.csv("hogehoge.csv",fileEncoding = "UTF-16",skip=1,sep="\t",header=T)
空白を空白として認識させておく
データを一部消したので公開した部分にはありませんが、何も書いていない列があったりします。そういう所は、欠損していると認識させておくほうが後で便利な場合があります。
コード
hogehoge <- read.csv("hogehoge.csv",fileEncoding = "UTF-16",skip=1,sep="\t",header=T,na.strings =c("","NULL")) heead(hogehoge)
結果の前の方*4
日 広告の状態 広告 広告文.1 広告文.2 表示URL リンク先URL 1 2016-07-13 アクティブ 広告名: 320x100_hogehoge.png; 320 x 100 <NA> <NA> URL <NA> 2 2016-07-06 アクティブ 広告名: 468x60_hogehoge.png; 468 x 60 <NA> <NA> URL <NA>
もっとちゃんとした情報は
その他にも、文字列が「factro」として認識されて「character」と認識されずに困ることがあります。
自分も苦労した経験があります。もっと早く知っておけばよかったと思います。
今回は必要としなかったその他の読み込み時のトラップはまだまだあります。読み込み時のオプションの指定方法についてはこちらのサイトがよくまとまっていました。
【備忘録】Python(x,y)からAnacondaに変えてみた
自分用のメモ。
Pythonはプロキシの関係でpipが使えないので、一通り最初から入っているPython(x,y)を使っていました。
しかし最近、Anacondaが話題なのと、Jupyterも使ってみたかったので、入れてみました。
これは、次に自分が入れるべき時に注意する点を整理したものです。
Anacondaをダウンロード
本家からWin版をダウンロード
Download Anaconda now! | Continuum
※参考にしたサイト
インストーラーでインストール
トラップは殆どなし。1箇所だけ、インストールするフォルダを自分で選ぶ場合に「もうその名前のフォルダはあるよ！」って怒られるので、フォルダを新しく作った場合は、マイコンピュータからアナログにフォルダを消します。
Jupyterの使い方
Jupyterをコマンドプロンプトから起動する場合、先ず、作業したいディレクトリにcdで移動しておき
jupyter notebook
で起動させるとブラウザが立ち上がる。
「New」⇒「Python2」でNotebookが開く。
実行して改行する為には、「Shift+Enter」で進む。
また、グラフを表示させる為には、
%matplotlib inline
と打っておく必要がある。
参考にしたサイト
Seabornを使ってみる
PythonにおけるRのggplot2的な存在の一つであるSeabornを使ってみた。
例のごとくcondaもpipも使えないので、直接ファイルをダウンロードする。
また、Anacondaの場合モジュールは
C:\(インストールしたフォルダ)\Lib\site-packages\
の中にあるので、「seaborn」など自分でフォルダを作りそこに解凍。
コマンドプロンプトを起動し、seabornを解凍したフォルダへ移動し、
python setup.py install
でインストール。
Seabornの使い方
seabortnを使うときは
import matplotlib.pyplot as plt
import seaborn as sns
と適当に名前をつけてインポートするだけでいい。
こんな感じ
参考にしたサイト
PDFにするにはpandocが必要
jupyterのダウンロード形式にPDFが選べるけれど、それにはpandocが必要らしい
公式サイトからインストーラーで導入する方が、Win機の場合はよいようです。
そもそものpandocの使い方
先ず、pandocの使い方を知る。
Inner Journeys: PandocとLuaLaTeXを使ったPDF出力でコードブロックをきれいに表示する
の通りに従ってみる。
1.メモ帳、Emacs等で/mdのMarkdown形式のファイルをUTF-8で作成・保存する。
#こんにちはPandoc
2.コマンドプロンプトで当該ファイルがある場所まで移動し、ファイル名を指定し、アウトプットのファイル名と使うTexのエンジンを指定する。日本語入りファイルの場合は、LuaLaTexが良い。*1
pandoc -f markdown input.md -s -o output.pdf --latex-engine=lualatex -V documentclass=ltjarticle
Pandocを使って、jupyterのファイルをTeXにする
本当は、jupyterの画面からdownload as でPDFにしたいのだけど、上手く行かなかったので、こちらのサイトを参考に、一度TeXにして、後からPDFに変換。
日本語が入っていないファイルの場合
jupyterの作業ファイル（今回の場合は、Untitled.ipynbという名前）があるフォルダまでコマンドプロンプトでカレントディレクトリを移動し、以下のコマンドを入力。
jupyter nbconvert --to latex Untitled.ipynb
すると、「Untitled.tex」というTeXファイルが出来るので、TeXのエディターで開いて後はいつも通りPDFにすればいい。((余談ですがこういう時、eclipseのpydevでは面倒。。。)
日本語が入っていない場合は、pdfLaTeXでも構わないが、日本語はエラーになってしまうため、先ほどの参考サイトに従ってテンプレートファイルをLuaLaTeXに変更して、新しくテンプレートファイルを作る。
今回は、ltjsarticle.tplxとして保存。
PATHの通った場所にテンプレートを置くか、作業フォルダに常にテンプレートをコピペするのかは悩ましい所。
jupyter nbconvert --to latex --template ltjsarticle.tplx Untitled2.ipynb
とコマンドプロンプトで打つとlatexファイルができるので、TeXのエディターでPDFにしてしまえばいい。
*1:普段はXeLaTeXの方が早くて好きなのだけど。
