|
Hi Seppe,
cannot get exactly the same (assuming normal distribution)
confint.cv <- function(x,alpha=.05, method="modmckay"){
# Calculate the confidence interval of the cv of the vector x
# Author: Kevin Wright
# See: Vangel, Mark. Confidence Intervals for a Normal Coefficient
# of Variation. American Statistician, Vol 15, No1, p. 21--26.
# x <- c(326,302,307,299,329)
# confint.cv(x,.05,"modmckay")
x <- na.omit(x)
v <- length(x)-1
mu <- mean(x)
sigma <- sqrt(var(x))
k <- sigma/mu
# CV > .33 may give poor results, so warn the user
if(k>.33) warning("Confidence interval may be very approximate.")
method <- casefold(method) # In case we see "McKay"
if(method=="mckay"){
# McKay method. See equation 15.
t1 <- qchisq(1-alpha/2,v)/v
t2 <- qchisq(alpha/2,v)/v
u1 <- v*t1
u2 <- v*t2
lower <- k/sqrt(( u1/(v+1) -1)*k*k + u1/v)
upper <- k/sqrt(( u2/(v+1) -1)*k*k + u2/v)
} else if (method=="naive"){
# Naive method. See equation 17.
t1 <- qchisq(1-alpha/2,v)/v
t2 <- qchisq(alpha/2,v)/v
lower <- k/sqrt(t1)
upper <- k/sqrt(t2)
} else {
# Modified McKay method. See equation 16.
method="modmckay"
u1 <- qchisq(1-alpha/2,v)
u2 <- qchisq(alpha/2,v)
lower <- k/sqrt(( (u1+2)/(v+1) -1)*k*k + u1/v)
upper <- k/sqrt(( (u2+2)/(v+1) -1)*k*k + u2/v)
}
ci <- cbind.data.frame(method = method, lower=lower,upper=upper, CV= k, alpha=alpha)
return(ci)
}
x=c(3.8,3.9,3.9,3.8,4.1,4.0,3.4,3.6,3.8,3.7)
results <- rbind(confint.cv(x), confint.cv(x, method = "mckay"), confint.cv(x, method = "naive"))
print(results)
gives
method lower upper CV alpha
1 modmckay 0.03617573 0.09632067 0.05263158 0.05
2 mckay 0.03618047 0.09641016 0.05263158 0.05
3 naive 0.03620185 0.09608475 0.05263158 0.05 —
Kind regards,
Mittyri |
Log in
Register
Ing. Helmut Schütz