ベイズの基礎：課題17-2 β分布の描画, polygon関数

課題17-1 β分布を事前確率に からの続き。


Activity 17-2: Employment Rate for Women

課題：米国で女性が働いている割合を求めたい。1988年の調査では「非就労率は5.5%」であった。

現在の就労率は1988年と同じレベルと考えられる。つまり、95% 近辺が妥当な率と想定する。しかし、この率は男女が含まれたものなので、女性だけの場合 95% よりも低いと考えられる。就労せずに家事・育児を選択することもあるから。


そんな女性の就労率 p を beta(18, 2) で表せると仮定すると、次のベータ分布となる：

curve(dbeta(x, 18,2), from=0, to=1, n=101, xlab="p", ylab="")
text(0.5,6,"beta(18,2) curve",cex=2)

　「連続型」なので plot 関数より、curve 関数ですよね (^^)b

以下は (c) から (h) の課題を解いたもの。

> pbeta(0.8,18,2) # pが80%以下の確率
[1] 0.08286623
> pbeta(0.9,18,2) # pが95%以下の確率
[1] 0.420265
> pbeta(0.9,18,2) - pbeta(0.8,18,2) # pが80%から90%の確率
[1] 0.3373987
> pbeta(0.85,18,2) - pbeta(0.77,18,2) # pが77%から80%の確率
[1] 0.151955
> 1 - pbeta(0.9,18,2) # pが90%以上の確率
[1] 0.579735
> 1 - pbeta(0.75,18,2) # pが75%以上の確率
[1] 0.9689926

pbeta 関数は「**%以下」、分布では原点から x 軸のプラス方向へ下位区間。「**%以上」の上位区間は  1 - pbeta(p) となる。


連続型分布に必須の描画

以降はテキストにはないが、連続型確率分布では欠かせない「領域塗りつぶし」を行った。この塗りつぶされた面積（積分値）が「確率」なのは「確率密度を例える」に記した。

beta(18, 2) の「90%信頼区間」をグラフで示す。

> curve(dbeta(x, 18,2), from=0, to=1, n=101, xlab="p", ylab="")

> qbeta(c(0.05,0.95),18,2)

> p <- qbeta(c(0.05,0.95),18,2)
> p
[1] 0.7736258 0.9809669
> x <- seq(p[1], p[2], by=(p[2]-p[1])/99)

> xval <- c(x,rev(x))

> yval <- c(rep(0,100),dbeta(rev(x),18,2))

> polygon(xval, yval, col=c("blue"), density=c(20), angle=c(45))
> text(0.7,6,"90% interval", col="blue")

両脇の区間を「赤の逆斜め線」で塗りつぶす。

> x1 <-seq(0,p[1],by=p[1]/49)

> x2 <-seq(p[2],1,by=(1-p[2])/49)

# 下位5%区間

> xval1 <- c(x1,rev(x1))

> yval1 <- c(rep(0,50),dbeta(rev(x1),18,2))

> polygon(xval1, yval1, col="red", density=20, angle=-45)

# 上位95%区間

> xval2 <- c(x2,rev(x2))

> yval2 <- c(rep(0,50),dbeta(rev(x2),18,2))

> polygon(xval2, yval2, col="red", density=20, angle=-45)

> text(0.7,5,"Lower 5%", col="red")

> text(0.7,4.5,"Upper 95%", col="red")

polygon関数

エリアを描画（塗りつぶす）方法で、polygon関数を使った。この関数がやっているのは、「描いた多角形の内面を脚色する」こと。

三角形で説明する。

> xval <- c(0,1,2,2,1,0)
> yval <- c(0,0,0,0,2,0)
> plot(xval, yval, col="red", lwd=10)

plot の引数は (x, y) 座標。

先の「赤点」を通る直線で描かれる「多角形」の内面を脚色：

> polygon(xval, yval, col=c("blue"), density=c(20), angle=c(45))

多角形の座標の作成は、rev, rep 関数が便利：

> x <- c(0,1,2)

> c(x,rev(x))

[1] 0 1 2 2 1 0

> c(rep(0,3),0,2,0)

[1] 0 0 0 0 2 0

課題17-3 β分布の決定方法 に続く。