希少サンプル

DNA鑑定試料がある。 試料の量が少ない。 1回の実験に全部使うのがよいのか、再実験用に取り分けておくのがよいのか?

基礎編

どうして再実験用に取り分けておきたいと感じるのか?

1回目で失敗しそう…

2回目は失敗しないのか?

失敗を想定して再実験する、とは

何回かに分けて実験する。

そのうちの1回でも成功すれば(成功して正しい結果が得られれば)、よし、とすることにする。

毎回、同じ確率で失敗するなら:

p.failure <- 0.5

plot(0:10,p.failure^(0:10),type="b",main="すべて失敗する確率",xlab="試行回数",ylab="確率")

複数回やれば、全部失敗する確率は、どんどん小さくなる。

\[ p^k \]

どうして分割をためらうのか?

全量を1回に使っても、ちゃんとした結果が出ないかもしれないのに、分割したら、全滅するかもしれない。 それよりは、できるだけ好条件の実験を1回だけやる方がよいのでは?

失敗の確率が1回の実験量と、たとえば次のような関係があるとする

\[ p.failure \times volume.exp = C \]

n.part <- 1:10
p1 <- 0.05
p.failure <- p1*n.part
plot(n.part,p.failure,xlab="分割数",ylab="確率",main="失敗する確率")

plot(n.part,p.failure^n.part,type="b",xlab="分割数",ylab="確率",main="すべて失敗する確率")

分割に従ってもっとどんどん成功率が下がる(失敗率が上がる)と…

n.part <- 1:10
p1 <- 0.01
p.failure <- p1*n.part^1.9
par(mfcol=c(1,2))
plot(n.part,p.failure,xlab="分割数",ylab="確率",main="失敗する確率")
plot(n.part,p.failure^n.part,type="b",xlab="分割数",ylab="確率",main="すべて失敗する確率")

par(mfcol=c(1,1))

理論編

分割してもしなくても同じ場合

容量 vの減少により、指数関数的に失敗率が下がるとき

v=0の失敗率が1

\[ f(v) = e^{-av} \]

v <- seq(from=0,to=3,length=100)
a <- 1
plot(v,exp(-a*v),type="l",xlab="容量",ylab="確率",main="失敗確率")

総量 v=1

分割してみる

(0.5,0.5) (0.2,0.8) (0.5,0.3,0.2)

# utility functions
my.g <- function(x,as,ps=0:(length(as)-1),scale=1){
    ret <- rep(0,length(x))
    for(i in 1:length(as)){
        ret <- ret + (as[i]/scale)*x^ps[i]
    }
    ret
}
my.div.g <- function(vs,as,ps=0:(length(as)-1),scale=1,log=FALSE){
    ret <- 0
    for(i in 1:length(vs)){
        ret <- ret + my.g(vs[i],as,ps,scale)
    }
    if(!log){
        ret <- exp(ret)
    }
    return(ret)
}

as <- c(-1)
ps <- c(1)
v <- seq(from=0,to=2,length=101)
scale <- 1
#plot(v,exp(my.g(v,as,ps,scale)),type="l")

2分割なら2回とも失敗する確率、 3分割なら3回とも失敗する確率を計算してみる。

vs <- c(0.5,0.5)
my.div.g(vs,as,ps,scale)
## [1] 0.3678794
vs <- c(0.8,0.2)
my.div.g(vs,as,ps,scale)
## [1] 0.3678794
vs <- c(0.5,0.3,0.2)
my.div.g(vs,as,ps,scale)
## [1] 0.3678794

分割の仕方によらず、すべて失敗する確率は同じ。

なぜなのか?

総量Vの試料をk回に分割する

\[ \sum_{i=1}^k v_i = V \]

各回の失敗率

\[ f(v) = e^{-av} \]

全回、失敗する確率は、以下の通り、分割の仕方によらずVのみによる

\[ \prod_{i=1}^k f(v) = \prod_{i=1}^k e^{-av} = e^{\sum_{i=1}^k v_i} = e^V \] \[ f(v) = e^{-av} \] という特別な失敗率と容量との関係があるとき、その分割の仕方によらず、全回、失敗する確率は総量によって決まる。

「指数関数」は「掛け算が足し算になる」ので、

「指数関数的に失敗率が決まるとき、分割の仕方によらず、全て失敗する確率は、使用する試料の総量によって決まる」

容量非依存性にも失敗するとすると

容量0のときの成功率0、容量が無限大のときの失敗率、K/(1+K)

\[ f(v) = \frac{e^{-av}+K}{1 + K} \]

K <- 1
plot(v,(exp(-v)+K)/(1+K),type="l",xlab="容量",ylab="確率",main="失敗確率")

二等分、三等分、四等分で、「全て失敗する確率」は変わってくる

vs <- rep(1/2,2)
prod((exp(-vs)+K)/(1+K))
## [1] 0.6452352
vs <- rep(1/3,3)
prod((exp(-vs)+K)/(1+K))
## [1] 0.6322156
vs <- rep(1/4,4)
prod((exp(-vs)+K)/(1+K))
## [1] 0.6257333

n等分するとき、nと「全て失敗する確率」との関係

n.v <- 1:10
out <- rep(0,length(n.v))
for(i in 1:length(n.v)){
  vs <- rep(1/n.v[i],n.v[i])
  out[i] <- prod((exp(-vs)+K)/(1+K))
}
plot(n.v,out,xlab="等分割数",ylab="確率",main="すべて失敗する確率")

失敗率関数を一般化してみる

指数関数がなかなか良い仕事をするので、指数関数を使う

ただし、失敗率関数を一般化したいので、指数を試料量の関数(特に多項式)で表すことにする

\[ f(v) = e^{g(v)} \]

\[ g(v) = -av \]

\[ g(v) = \sum_{i=1}^k a_i v^{p_i} \]

g(v) が単純な場合

\[ g(v) = a \times v^p \]

p=0の場合

\[ g(v) = a \times v^0 = a \]

as <- c(-0.5)
ps <- c(0)
v <- seq(from=0,to=2,length=101)
scale <- 1
plot(v,exp(my.g(v,as,ps,scale)),type="l",xlab="容量",ylab="確率",main="失敗確率")

等分割するとしたら

n.v <- 1:10
out <- rep(0,length(n.v))
for(i in 1:length(n.v)){
  vs <- rep(1/n.v[i],n.v[i])
  out[i] <- my.div.g(vs,as,ps,scale)
}
plot(n.v,out,xlab="分割数",ylab="確率",main="全て失敗する確率")

p=1の場合

\[ g(v) = a \times v^1 = av \]

as <- c(-1)
ps <- c(1)
v <- seq(from=0,to=2,length=101)
scale <- 1
plot(v,exp(my.g(v,as,ps,scale)),type="l",xlab="容量",ylab="確率",main="失敗確率")

等分割する

これは、「失敗率関数が単純な指数関数の場合」=「全部失敗する確率は分割の仕方によらない」

n.v <- 1:10
out <- rep(0,length(n.v))
for(i in 1:length(n.v)){
  vs <- rep(1/n.v[i],n.v[i])
  out[i] <- my.div.g(vs,as,ps,scale)
}
plot(n.v,out,xlab="分割数",ylab="確率",main="少なくとも1回は成功する確率")

p=2の場合

\[ g(v) = a \times v^2 = av^2 \]

as <- c(-1)
ps <- c(2)
v <- seq(from=0,to=2,length=101)
scale <- 1
plot(v,exp(my.g(v,as,ps,scale)),type="l",xlab="容量",ylab="確率",main="失敗確率")

等分割する

n.v <- 1:10
out <- rep(0,length(n.v))
for(i in 1:length(n.v)){
  vs <- rep(1/n.v[i],n.v[i])
  out[i] <- my.div.g(vs,as,ps,scale)
}
plot(n.v,out,xlab="分割数",ylab="確率",main="すべて失敗する確率")

p=0.5の場合

\[ g(v) = a \times v^{0.5} \]

as <- c(-1)
ps <- c(0.5)
v <- seq(from=0,to=2,length=101)
scale <- 1
plot(v,exp(my.g(v,as,ps,scale)),type="l",xlab="容量",ylab="確率",main="失敗確率")

等分割する

n.v <- 1:10
out <- rep(0,length(n.v))
for(i in 1:length(n.v)){
  vs <- rep(1/n.v[i],n.v[i])
  out[i] <- my.div.g(vs,as,ps,scale)
}
plot(n.v,out,xlab="分割数",ylab="確率",main="すべて失敗する確率")

1,2,3分割

分割数を1,2,3のいずれかにする。

分割はいろいろな割合にする。

\(g(v) = a\times v^2\) の場合

\[ g(v) = a \times v^2 = av \]

分割する方がよく、分割数が多い方がよく、均等分割が一番よい。

library(MCMCpack)
## Loading required package: coda
## Loading required package: MASS
## ##
## ## Markov Chain Monte Carlo Package (MCMCpack)
## ## Copyright (C) 2003-2017 Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park
## ##
## ## Support provided by the U.S. National Science Foundation
## ## (Grants SES-0350646 and SES-0350613)
## ##
library(mwaytable)
n.pt <- 10^4
Vs <- rdirichlet(n.pt,rep(1,3))
cv <- CategoryVector(3)
xy <- Vs %*% cv
as <- c(-1)
ps <- c(2)
v <- seq(from=0,to=2,length=101)
scale <- 1

out <- rep(0,n.pt)
for(i in 1:n.pt){
  vs <- Vs[i,]
  out[i] <- my.div.g(vs,as,ps,scale,log=TRUE)
}


out.st <- (out-min(out))/(max(out)-min(out))
plot(xy,pch=20,col=rgb(out.st,1-out.st,1),xlim = c(-1,1),ylim=c(-1,1))

library(rgl)
## Warning: package 'rgl' was built under R version 3.2.5
library(knitr)
knit_hooks$set(rgl = hook_rgl)
plot3d(xy[,1],xy[,2],out)