在看binom.test()的时候,看到了confidence interval的计算。
看到源码发现了关于置信区间的计算,利用到了二线分布和贝塔分布的关系。
## Determine p s.t. Prob(B(n,p) >= x) = alpha.
## Use that for x > 0,
## Prob(B(n,p) >= x) = pbeta(p, x, n - x + 1).
p.L <- function(x, alpha) {
if(x == 0) # No solution
0
else
qbeta(alpha, x, n - x + 1)
}
## Determine p s.t. Prob(B(n,p) <= x) = alpha.
## Use that for x < n,
## Prob(B(n,p) <= x) = 1 - pbeta(p, x + 1, n - x).
p.U <- function(x, alpha) {
if(x == n) # No solution
1
else
qbeta(1 - alpha, x + 1, n - x)
}
#置信区间
CINT <- switch(alternative,
less = c(0, p.U(x, 1 - conf.level)),
greater = c(p.L(x, 1 - conf.level), 1),
two.sided = {
alpha <- (1 - conf.level) / 2
c(p.L(x, alpha), p.U(x, alpha))
})
对于Prob(B(n,p) >= x) = pbeta(p, x, n - x + 1).和Prob(B(n,p) <= x) = 1 - pbeta(p, x + 1, n - x).的推到过程不了解,希望有懂的朋友指点一下 :-)