R/binomci.R

binomci <-
function(x=sum(y),nn=length(y),y=NULL,n=NA,alpha=.05){
#  Compute a 1-alpha confidence interval for p, the probability of
#  success for a binomial distribution, using Pratt's method
#
#  y is a vector of 1s and 0s.
#  x is the number of successes observed among n trials
#
if(!is.null(y)){
y=elimna(y)
nn=length(y)
}
if(nn==1)stop("Something is wrong: number of observations is only 1")
n<-nn
if(x!=n && x!=0){
z<-qnorm(1-alpha/2)
A<-((x+1)/(n-x))^2
B<-81*(x+1)*(n-x)-9*n-8
C<-(0-3)*z*sqrt(9*(x+1)*(n-x)*(9*n+5-z^2)+n+1)
D<-81*(x+1)^2-9*(x+1)*(2+z^2)+1
E<-1+A*((B+C)/D)^3
upper<-1/E
A<-(x/(n-x-1))^2
B<-81*x*(n-x-1)-9*n-8
C<-3*z*sqrt(9*x*(n-x-1)*(9*n+5-z^2)+n+1)
D<-81*x^2-9*x*(2+z^2)+1
E<-1+A*((B+C)/D)^3
lower<-1/E
}
if(x==0){
lower<-0
upper<-1-alpha^(1/n)
}
if(x==1){
upper<-1-(alpha/2)^(1/n)
lower<-1-(1-alpha/2)^(1/n)
}
if(x==n-1){
lower<-(alpha/2)^(1/n)
upper<-(1-alpha/2)^(1/n)
}
if(x==n){
lower<-alpha^(1/n)
upper<-1
}
phat<-x/n
list(phat=phat,ci=c(lower,upper),n=n)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.