読者です 読者をやめる 読者になる 読者になる

tak0kadaの何でもノート

発声練習、生存確認用。

医学関連は 医学ノート

2項分布、ベータ分布をRとpythonで描いてみる

R 統計 python matplotlib scipy.stats 分布

2項分布

  • R $$ P(X=k) = \begin{pmatrix} n\\ k\\ \end{pmatrix} p^{k}(1-p)^{n-k} $$
library("rgl")

#サイコロを投げる最大の回数
N <- 200
#成功率
p <- 0.3

#投げる回数をx、成功する割合をy、その場合の確率をzとする
x <- y <- z <- rep(0, 1/2*(N+1)*(N+2)-1)
head <- 1
for (n in 1:N){
  x[head:(head+n)] <- n
  y[head:(head+n)] <- (0:n) / n
  z[head:(head+n)] <- dbinom(0:n, n, p)
  head <- head + n + 1
}

plot3d(x, y, z, pch=20, col="blue", xlab="n_trial", ylab="success_rate", zlab="p")
rgl.snapshot("binom.png")

Nが大きくなるとそれぞれの仮説ごとの確率は小さくなるが0.3付近に収束して行っているのが分かる。

f:id:kuyata:20140617212251p:plain

from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolskits.mplot3d import Axes3D
from scipy.stats import binom

N = 200
p = 0.3

x = np.hstack([[n]*(n+1) for n in range(1, N)])
y = np.hstack([np.arange(n+1)/n for n in range(1, N)])
z = np.hstack([binom.pmf(np.arange(n+1), n, p) for n in range(1, N)])

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x, y, z, c="blue", marker=".", linewidrh="0")
ax.set_xlabel('n_trial')
ax.set_ylabel('success')
ax.set_zlabel('p')

plt.show()

f:id:kuyata:20140618125220j:plain

軸のラベルの設定忘れれたorz

ベータ分布

例えばサイコロを振って$m $回、$n$回表と裏が出た時、 $\alpha = m+1, \beta = n+1$として、

$$ p(x, \alpha, \beta) = \dfrac{1}{B(\alpha, \beta)}x^{\alpha-1}(1-x)^{\beta-1} $$ $$ B(\alpha, \beta) = \dfrac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} $$

wikipediaと同じ図を描いてみる

  • R
library("rgl")

#グラフは折れ線グラフで近似する。分割の個数。
n_dev = 100

alpha <- c(0.5, 5, 1, 2, 2)
beta <- c(0.5, 1, 3, 2, 5)
f <- function(a, b) dbeta(seq(from=0, to=1, length=n_dev), a, b)
z <- mapply(f, alpha, beta)

#グラフを重ねるのにpar(new=T)を使う
for (i in 1:5){
  plot(seq(from=0, to=1, length=n_dev) , z[,i], xlim=c(0,1), ylim=c(0,3), col=i, type='l')
  par(new=T)
}

Rはimport文を書かなくても結構使えるのでさくっと計算するときにも便利である。

f:id:kuyata:20140617011108p:plain

from scipy.stats import beta
import numpy as np
import matplotlib.pyplot as plt

#n_devは100だとa=b=0.5のときにx=0,1付近で無限大になるのが余り分かりやすくなかった
n_dev = 1000
a = (0.5, 5, 1, 2, 2)
b = (0.5, 1, 3, 2, 5)

#beta.pdf(x, alpha, beta)を使う
x = np.linspace(0, 1, num=n_dev)
f = lambda a, b: beta.pdf(x, a, b)
y = np.array([beta.pdf(x, a, b) for a, b in zip(a, b)])

for i in range(5):
    plt.plot(x, y[i])

plt.show()

f:id:kuyata:20140617015019p:plain

pdfなのでこれ単体で意味を持つわけではないが、例えば尤度を計算するときにはpmfをビンごとに積分して求める、あるいは一次近似で良いならpdfをそのまま使って後で正規化することもできる。