• きぬいと
  • TokyoRで毎月発表してたような人。
  • 7分枠で20分しゃべる人
• もともと仙台で6年間勉強とか研究とかしてた
  • ここのために仙台に来た
  • 終わったら帰ります
• 今は調査屋
  • 主なお仕事: 市場調査とマーケティングとやらかしのリカバリ
  • ウチで働きたい人募集してます。

• 明日以降Twitterで公開します
  • Rmarkdownでつくってます
• 内容は
  • 7月のTokyoRで「話せなかった」お話です。

• 統計学寄りのお話
  • 確率「分布」のお話
    • 「確率分布に従う」ことについてはこちら
• 「正規分布に従う」ってよく聞く
• しかし我々は「正規分布」ってよくわからない
• なんか「確率分布の祖」みたいなところに行きたい
  • 人生で必ず一回は目にするであろう統計屋の地図があります。

• by: http://www.math.wm.edu/~leemis/chart/UDR/UDR.html

• 言われてみれば代表値なんもわからんくないですか？
  • せめて平均だけでも持って帰ってほしい
  • あとついでに分散も

• 観測値の総和を試行回数で割る。
  • 今回はコイントス
    • 表が出たら1、裏なら0と記録する係を雇う
• きぬいと: 人望がないのでRに頼る
  • コイントス100回の試行を100回繰り返す。
  • すでに「分布」に頼ってしまっている(run-ifではなくr-unifなんだね)

X <- numeric(100)
for(try in 1:100){
  set.seed(20+try) # seed値を毎回変える
  X[try] <- sum(round(runif(100, min = 0, max = 1)))
}
mean_X <- sum(X)/length(X)
mean_X

## [1] 49.76

• つまりコイントス100回で49.76回表が出そうだ、ということになる。
• 「歪みのないコインでコイントスをしたら半分くらいは表になるんじゃね？」は合ってそう……
• 平均は強い。
  • 数多の機械学習モデルが複雑な計算を経て出した結果でも、平均に勝てないこともある。
  • AIは大体何かしらの条件付き平均だってTokyoR設立者の人も言ってた

• なんだかんだ「ズレ」る
  • どうにかズレの度合いの情報も知りたい。
  • 平均からの差がどれだけ離れているかを見ればよいのでは？

• 覚悟するための数字実現した値の、平均からの離れ具合
• \(Var(X) = \frac{1}{N}\sum^{N}_{i}(X_i - \overline{X})^2\)
  • 2乗する理由: 足し上げる時にマイナスになると悲しいから
• 正の2乗根をとると標準偏差(単位が平均と同じ)になる。

Var_X <- sum((X-mean_X)^2)/length(X)
sqrt(Var_X)

## [1] 5.05197

• 大体5回ぐらい散らばりそう。
  • 45〜55回の間で答えられればいい……？
• もしかしたら45に寄ってるかも……あるいはその逆かも……
• どれくらい直感から「歪んで」いるんだろうなぁ……？

• ワイド？
  • 散らばり方のバランスをみる。
  • \(\alpha = \frac{1}{N}\sum^{i}_{N}(\frac{X_i-\mu}{var(X)})^3\)
  • 分子が正負どっちもとりうる
  • 平均を真ん中においたとき、左右どっちに「歪み」がちか？
• みんな辞書登録しとこうね(ゆがみどで変換)

Skew_X <- sum(((X-mean_X)/Var_X)^3)/length(X)
Skew_X

## [1] 0.001586152

• 歪みねえな

• これらはすべて“モーメント”という概念で定式化できるらしい
  • 詳しくは未来のきぬいとに託す
• データの歪みとかは具体的に数値としてみることは多くない印象
  • ヒストグラムを見て判断することが多い
  • でもこの世界はまだ「データの可視化」以前の世界なので……

• 対象者を増やしてコイントスをしてもらったらどうなるんだろう

X <- numeric(1000)
for(try in 1:length(X)){
  set.seed(20+try) # seed値を毎回変える
  X[try] <- sum(round(runif(100, min = 0, max = 1)))
}
mean_X <- sum(X)/length(X)
mean_X

## [1] 49.915

X <- numeric(10000)
for(try in 1:length(X)){
  set.seed(20+try) # seed値を毎回変える
  X[try] <- sum(round(runif(100, min = 0, max = 1)))
}
mean_X <- sum(X)/length(X)
mean_X

## [1] 50.0091

• 平均は大体試行回数の半分っぽい
• 対象者数に依存せず、常に一定の「割合」で表が出てきそう
  • この割合をパラメータ(確率)として設定しちゃおう(\(p\))
  • さらに、事前にわかる試行回数(\(n\))も使っちゃおう
  • 表の出る回数はやってみないとわからないが、\(n\)を越えることはないとか、そういうのはわかるなあ
• つまり、表の出る回数に関する情報を、\(n\)と\(p\)を使ってうまいこと計算できたらいいね

• 表が出るときの「出方のパターン」は何通りか？
  • 表→表→裏とか裏→表→表とか。いろいろある
  • みんな勉強したはずの「組合せ」を思い出す
• 50回表が出るとき、コイントスをした順に並べた際の並べ方のパターンは…… -　\({}_{100} \mathrm{C}_{50}\)通り……
  • 数字的には\(1.0089134\times 10^{29}\)通り……こわ……何……それ……

• 話を簡単にするために最初の半分は連続で表が出たとしましょう。
• コイントスは練習してもうまくならないことにすると(独立)
• 一定の割合\(p = 0.5\)で表が出るのであれば、裏が出る割合は\(1-p = q = 0.5\)
  • 最初の50回表が出て、残りで裏が出る確率は\(p^{50}q^{50}\)
  • 「表が50回出て、裏が50回出る」パターンは、全部で\(1.0089134\times 10^{29}\)
  • このパターン数だけ\(p^{50}q^{50}\)を足すことで「50回表の出る確率」が決まる

• \(P(表回数=50)={}_{100} \mathrm{C}_{50}p^{50}(1-p)^{50}=0.0795892\)
• えっ、ちっさ
  • 故に「45〜55回出るんじゃない？」的な「範囲」で確率を定義してあげる
  • そうすると大体\(0.728747\)

• コイントス\(n\)回をするとき、実現する回数\(k\)の分布
  • \(P(k) = {}_{n} \mathrm{C}_{k}p^{k}(1-p)^{n-k}\)
• これが二項分布

• \(n\)回中、表の出る回数\(k\)全パターン(0〜\(n\))それぞれに、生起確率をかける
  • \(\mu = \sum^{N}_{i}k_iP(k_i)\)
• 二項分布によって運命づけられたコイントスの回数の平均が出てくる

BIN <- function(n,k,p){choose(n,k)*(p^k)*((1-p)^(n-k))} #確率分布の定義まんま
n_ <- 100; k_ <- c(0:100); p_ <- 0.5 # 確率0.5くらいでみてみようね
mu <- sum(k_*BIN(n_,k_,p_)); print(mu)

## [1] 50

• \(\mu = np\)なんですよ。

• 平均が出れば分散は簡単
• \((回数-平均)\)の生起確率も回数の生起確率に依存します。
  • \(var(k) = \sum^{N}_{i}(k_i-\mu)^2P(k_i)\)

diff_k <- (k_-mu)^2
Var_X <- sum(diff_k*BIN(n_,k_,p_)); print(Var_X);sqrt(Var_X)

## [1] 25

## [1] 5

• \(25 = 100*0.5*0.5 = np^2\)？
• 実際は\(np(1-p)\)(標準偏差は\(\sqrt{np(1-p)}\))

• 二項検定で使う。
  • 自社商品AとBの売上にちゃんと差があるのか？
• ロジスティック回帰で使う
  • 二項分布へのリンク関数がロジスティック関数という解釈
    • glm(y ~ ., family = binomial("logit"))←ほら
  • プロビット関数でも可
  • (機械学習的な話なら)シグモイドでもいいよね。

• 次なる知的好奇心が生まれる -　\(n\)とか\(p\)とか、動かしてみたくないですか？
  • 本当におっきい\(n\)とかでも、小さい\(p\)とかで振る舞い方は変わらないんだろうか？
• もっと一般化したい
  • 回数以外に使えないのは厳しい。いろいろな値をとることも考えられる

• 続きはまた、そのうち。