二項分布・ベータ分布、ポアソン分布・ガンマ分布

作者: ryamada
発売日: 2014/08/03
メディア: Kindle版
この商品を含むブログ (1件) を見る
事後分布推定、プライアは、とか、色々あるけれど、ひとまず、事前情報なしでデータがあったときに、Rを叩いてみることだけを目指しての資料
Rmdファイルは以下→勝手にhtml化・epub化するにはこちらを参照
htmlファイルはこちらにゲストログインで"16 September - 22 September"のテキストから。
# 事後分布推定 Estimation of posterior density function

## 成功確率、尤度 Success rate, likelihood

> 希少難病Xの患者を治療法Aで、５人治療したことがあります。３人で成功、２人で失敗でした。治療法Ａの成功率 p.a.sample はどれくらい？

> Assume a rare intractable disease, X, and you have treated five patients of X with a treatment A; 3 successes and 2 failures. How high is A's success rate p.a.sample ?

```{r}
# 成功回数、失敗回数 No. successes and failures
n.s <- 3
n.f <- 2

# 「成功した」割合 Fraction of successes among trials
p.a.sample <- n.s/(n.s + n.f)
p.a.sample
```

> ６人目の患者さんの治療成功率はいくつ？

> How likely would your 6-th patient be successfully treated?

こんな風に考えてみる。Think this way:

> Aの成功率が、たとえば、pa = 0.5 だとしたら、５人治療して３人成功する確率 pr0.5はいくつか？

> When A's success rate is pa = 0.5, what is the probability to have 3 successes among 5 trials, pr0.5?

まず準備。Preparation.

```{r}
# 3! = 1*2*3
factorial(3)
factorial(1:5)
# Combination 3 out of 5: 5!/(3!2!)
factorial(5)/(factorial(3)*factorial(2))
choose(5,3)
choose(n.s+n.f,n.s)
pa <- 0.5
pr.0.5 <- choose(n.s+n.f,n.s) * pa^n.s*(1-pa)^(n.f)
pr.0.5
```

> Aの成功率が、pa=0.7だとしたらどうなるか？

> Assume pa=0.7, then?

```{r}
pa <- 0.7
pr.0.7 <- choose(n.s+n.f,n.s) * pa^n.s*(1-pa)^(n.f)
pr.0.7
```

pa = 0.5, 0.7ではどちらが「ありそう」か。
pa = 0.5 or 0.7, which is more likely?

pr.0.5, pr.0.7 はこの観測に対する、pa=0.5, 0.7の尤度

pr.0.5, pr.0.7 are likelihood of pa=0.5, 0.7 with this observation, respectively.

```{r}
pr.0.5
pr.0.7
# ratio of likelihood: "likelihood ratio"
pr.0.5/pr.0.7
```

> Aの成功率pa は0.5かもしれないし0.7かもしれないし、それ以外の色々な値かもしれないので、尤度をそれらについて調べることにする
> pa can be 0.5 or 0.7 and other various values. Let calculate likelihood of the values.

```{r}
pa <- seq(from=0,to=1,length=100)
pr <- choose(n.s+n.f,n.s) * pa^n.s*(1-pa)^(n.f)
plot(pa,pr,type="l")
```

## 便利な分布、ベータ分布 Useful distribution, beta distribution

```{r}
pr.beta <- dbeta(pa,n.s+1,n.f+1)
plot(pa,pr.beta,type="l")
```

２つのコマンド、２つのグラフを比べてみよう。何が違って何が同じ？

## ベータ分布についての知識 Knowledge on beta distribution

今の疑問はこれ。
What you want to know is this.

> ６人目の患者さんの治療成功率はいくつ？

> How likely would your 6-th patient be successfully treated?

http://ja.wikipedia.org/wiki/ベータ分布

http://en.wikipedia.org/wiki/Beta_distribution

記事を読まずに情報カラムを読む。
Do not read sentences but read information column.

## 期待値・平均値、最頻値、中央値　Expected value/mean, mode and median

```{r}
alpha <- n.s+1
beta <- n.f+1
me <- alpha/(alpha+beta)
me
(n.s+1)/((n.s+1)+(n.f+1))
md <- (alpha-1)/((alpha-1)+(beta-1))
md
n.s/(n.s+n.f)
pr.beta <- dbeta(pa,n.s+1,n.f+1)
plot(pa,pr.beta,type="l")
abline(v=me,col=2)
abline(v=md,col=3)
```

中央値の式はちょっと変。The formula of mode is a bit wierd.

乱数使って求めてみる。 We can get the value through random number generation.

```{r}
n <- 100000
r <- rbeta(n,n.s+1,n.f+1)
quantile(r,0.5)
(alpha-1/3)/(alpha+beta-2/3)
```

## ２つの治療法 Two therepies

> 希少難病Xの患者を治療法Aで、５人治療したことがあります。３人で成功、２人で失敗でした。治療法Bの結果は、７人治療して４人成功でした。


> Assume a rare intractable disease, X, and you have treated five patients of X with a treatment A;3 successes and 2 failures. With another treatment B; 4 successes and 3 failures.

```{r}
n.s.a <- 3
n.f.a <- 2
n.s.b <- 4
n.f.b <- 3
```

２つの治療法の成功率の期待値はいくつか？ Answer expected success rate of two treatments.

```{r}

n.iter <- 10000
n.attempt <- 1000
pa <- rbeta(n.iter,n.s.a+1,n.f.a+1)
pb <- rbeta(n.iter,n.s.b+1,n.f.b+1)

(n.s.a+1)/(n.s.a+n.f.a+2)
(n.s.b+1)/(n.s.b+n.f.b+2)
```

> ２つの実験方法A,Bがある。予備実験で、Aは3回成功、2回失敗、Bは4回成功、3回失敗した。本実験をする。本実験は1000回まとめて実施する。1000回のうち400回以上、成功してくれないと、結果の評価ができない。Aを選ぶかBを選ぶか

> There are two experimental methods A and B. You performed preliminary experiments and with A, 3 successes and 2 failures, with B, 4 successes and 3 failure. You have to do 1000 exepriments all together. Among 1000, you have to have at least 400 successes so that you can step further. Which method would you select, A or B, or something else?


```{r}
# Strategy 1: A
res1 <- rep(0,n.iter)
for(i in 1:n.iter){
  s <- sample(0:1,n.attempt,prob=c(1-pa[i],pa[i]),replace=TRUE)
  res1[i] <- sum(s)
}
# Strategy 2: B
res2 <- rep(0,n.iter)
for(i in 1:n.iter){
  s <- sample(0:1,n.attempt,prob=c(1-pb[i],pb[i]),replace=TRUE)
  res2[i] <- sum(s)
}
# Strategy 3: A and B both
res3 <- rep(0,n.iter)
for(i in 1:n.iter){
  num.a <- n.attempt/2
  num.b <- n.attempt - num.a
  sa <- sample(0:1,num.a,prob=c(1-pa[i],pa[i]),replace=TRUE)
  sb <- sample(0:1,num.b,prob=c(1-pb[i],pb[i]),replace=TRUE)
  res3[i] <- sum(sa) + sum(sb)
}

# Strategy 4: Assign A or B at random
#res4 <- rep(0,n.iter)
#for(i in 1:n.iter){
#  a.or.b <- sample(c("a","b"),n.attempt,replace=TRUE)
#  num.a <- length(which(a.or.b=="a"))
#  num.b <- n.attempt - num.a
#  sa <- sample(0:1,num.a,prob=c(1-pa[i],pa[i]),replace=TRUE)
#  sb <- sample(0:1,num.b,prob=c(1-pb[i],pb[i]),replace=TRUE)
#  res4[i] <- sum(sa) + sum(sb)
#}


h <- hist(c(res1,res2,res3))
hist(res1,add=TRUE,col=2,density=9)
hist(res2,add=TRUE,col=3,density=11)
hist(res3,add=TRUE,col=4,density=13)
legend(0,3000,c("1","2","3"),text.col=2:4)
#hist(res4,add=TRUE,col=5,density=13)
matplot(apply(cbind(res1,res2,res3),2,sort),type="l",col=2:4,ylab="No. successes")
abline(h=400)
legend(0,1000,c("1","2","3"),text.col=2:4)
```

Q 今日の宿題。Today's assignment.

AとBとを同じ回数実施するという戦略もありだが、AとBとに異なる回数を割り当てることにしてもよい。たとえば、Aに(0,100,200,...,900,1000)回を割り当て、残りをBにしたとして、どの場合が400以上の成功を得る確率を最大にするのはどの場合になるかを上記のようなシミュレーションにて調べよ。可能ならば、もっと場合分けを増やし、場合分けの仕方と、400以上の成功を得る確率との関係をプロットにせよ。

As we did above, you can take a strategy to assign half to A and the rest to B. You can think of different assignment strategy. For example you can assign (0,100,200,...,900,1000) to A and the rest to B. You can check probability to get more than 400 successes for each strategy. Answer the best assingment strategy by simulation. If you can do more, check more values of assigment to A and make a plot that shows relation between the number of experiment A and probability with 400 or more successes.

## ポアソン分布とガンマ分布 Poisson distribution and Gamma distribution

使わない日も多い試薬(もしくは使い捨ての装置)Xがある。Xは受注生産性だと言う。発注から納品までd日かかる。
Xは保管コストが高いので、基本的には、ストックはしたくないが、いくつかはストックせざるをえない。
このようにすることにした。
> n個をストックしておき、必要になったら使う。使うたびに１つ、補充分を発注する。

Xが必要となったのはいつだったかという記録を基に、ストック数nの数を検討したい。

There is a drug (or disposable device) which is not frequently used, X. X is produced on demand. It takes d days from the oder to replacement.
Because of expensive cost to keep X at hospitals, you don't want to keep more than neccesary.
You have a plan as:
> Keep n Xs at hand. When you use one, you place an order of 1 X.

Based on the record of usage of X in the past, you want to decide n.

## ランダムに起きるまれなイベント Rare events

平均して1日にk回起きるとする。
In average, k events per day.
ポアソン乱数を使ってデータを作ってみる。
Let's generate a data set with random value generating function in poisson distribution.

```{r}
k <- 1.5
ev.1 <- rpois(365,k)
plot(1:365,ev.1,type="h")
table(ev.1)
hist(ev.1)
```

30日分のデータがこのように得られたとして、逆にどのようなポアソン分布であるかを推定してみる。

Let's estimage back the poisson distribution based on this generated data set.

ポアソン分布のパラメタkは0より大きい実数(１よりも大きくなるからベータ分布ではだめ)。
The parameter of poisson distribution is more than 0 and can be more than 1, meaning beta distribution is not good for this parameter.

色々な値で尤度を計算してみる。

Let's calculate likelihood for various k values.

```{r}
k <- 1.5
# 30 days
ev.2 <- rpois(30,k)
table(ev.2)

ks <- seq(from=0.1,to = 10,length=1000)
like <- rep(0,length(ks))
max.events <- max(ev.2)
tab <- rep(0,max.events+1)
for(i in 1:length(tab)){
  tab[i] <- length(which(ev.2==i-1))
}
for(i in 1:length(ks)){
  # 尤度の相対比だけを気にするときは、組合せ項であるchoose()の項は不要。なぜなら、すべての尤度の計算で共通だから
  #The term of choose() is common for all k values. Therefore you can omit its calculation when you compare likelihoods from various k values.
  like[i] <- 0
  for(j in 1:length(tab)){
    # 尤度は掛け算だけれど、対数を取って足しても同じ
    # When you calculate likelihood, you use multiplication, but you can sum up logs instead.
    like[i] <- like[i] + tab[j] * log(dpois((j-1),ks[i]))
  }
}
par(mfcol=c(1,2))
plot(ks,like)
abline(v=1.5)
plot(ks,exp(like))
abline(v=1.5)
par(mfcol=c(1,1))
```

成否を扱ったとき、その成功率はベータ分布で扱えた。
ポアソン分布に従うと考えるとき、観察値からポアソン分布のパラメタがどんな値かを推定するにはガンマ分布を使う。

When estimating success rate, we used beta distribution. For poisson distribution, we use gamma distribution.

```{r}
n.total.events <- sum(ev.2)
n.days <- length(ev.2)
d.gamma <- dgamma(ks,shape=n.total.events+1,scale=1/n.days)
plot(ks,d.gamma)
abline(v=1.5)
# 尤度計算結果との一致を見る
# Correspondence between likelihood calculated and gamma distribution
plot(exp(like),d.gamma)
```

元の問題設定に戻る。
Now we should be back to the original problem.

> ストック数をnとし、補充にd日かかるとして、5年間で在庫切れが起きるかをシミュレーションしてみる。

> When we stokc n and replacement takes d days, how many times do we find out of stock over five years?

```{r}
# In case of No. stock being 2
ns <- 20
# Days for replacement
d <- 7
# Years to follow
y <- 5
# Set poisson parameter with random generator of gamma distribution
r.gamma <- rgamma(1,shape=n.total.events+1,scale=1/n.days)

r <- rpois(365*y, r.gamma)
# replacements are d days after events
replace <- rep(0,length(r))
replace[(d+1):length(r)] <- r[1:(length(r)-d)]
#cumulative calculation tells you stock condition
cum.consumption <- cumsum(r)
cum.replace <- cumsum(replace)
  
over.consumption <- cum.replace - cum.consumption
plot(over.consumption,type="l")
abline(h=ns+0.1,col=2)
length(which(over.consumption <= -ns))
```

５年間に何度までなら在庫切れを許すか？

How many days would you allow out-of-stock? 
```{r}
# No. stocks are parameterize
ns <- 1:40
# Days for replacement
d <- 7
# Years to follow
y <- 5
# Parameters for poisson distribution based on 30 days record
n.iter <- 10000
r.gamma <- rgamma(n.iter,shape=n.total.events+1,scale=1/n.days)

# results
res <- matrix(0,n.iter,length(ns))

for(i in 1:n.iter){
  r <- rpois(365*y, r.gamma[i])
  # replacements are d days after events
  replace <- rep(0,length(r))
  replace[(d+1):length(r)] <- r[1:(length(r)-d)]
  #cumulative 
  cum.consumption <- cumsum(r)
  cum.replace <- cumsum(replace)
  
  over.consumption <- cum.replace - cum.consumption
  # Check for various stock numbers
  for(j in 1:length(ns)){
    res[i,j] <- length(which(over.consumption <= -ns[j]))
  }
}
boxplot(res)
```