rm(list=ls(all=TRUE)) #### R code for the 95% confidence intervals for SD and CV #### Alternative Confidence Interval for Variability #### Parameters in the Normal Distribution with #### Applications to Stock Exchange Index Data Set #### contact: patarawan.s@gmail.com ############################################### #### Use function for 95% CIs for SD and CV ############################################### estimation = function(x, n){ ### Estimated SD alpha = 0.05 xbar = mean(x) s2 = var(x) s = sd(x) var.s = s2*(1- (2/(n-1))*(gamma(n/2)/gamma((n-1)/2))^2) s.ub = s*sqrt(n-1)*gamma((n-1)/2)/(sqrt(2)*gamma(n/2)) var.s.ub = (((n-1)/2)* (gamma((n-1)/2)/gamma(n/2))^2 - 1)*s.ub^2 cat("S =", s, ", S.ub =",s.ub ,'\n') ### CI_trad for sigma using original method in the book l.sigma.ori = sqrt(n-1)*s/sqrt(qchisq(1-alpha/2, df=n-1 )) u.sigma.ori = sqrt(n-1)*s/sqrt(qchisq(alpha/2, df=n-1 )) length.sigma.ori = u.sigma.ori-l.sigma.ori cat("CItrad =", l.sigma.ori, u.sigma.ori, ", length =", length.sigma.ori ,'\n') ### CI_adj (aa=1) for sigma using adjusted original method l.sigma.adj = sqrt(n-1)*s.ub/sqrt(qchisq(1-alpha/2, df=n-1 )) u.sigma.adj = sqrt(n-1)*s.ub/sqrt(qchisq(alpha/2, df=n-1 )) length.sigma.adj = u.sigma.adj-l.sigma.adj cat("CIadj =", l.sigma.adj, u.sigma.adj, ", length =", length.sigma.adj ,'\n') ### CI_adj (aa=2) aa = 2 l.sigma.adj2 = sqrt(n-aa)*s.ub/sqrt(qchisq(1-alpha/2, df=n-1 )) u.sigma.adj2 = sqrt(n-aa)*s.ub/sqrt(qchisq(alpha/2, df=n-1 )) length.sigma.adj2 = u.sigma.adj2-l.sigma.adj2 cat("CIadj.aa=2 =", l.sigma.adj2, u.sigma.adj2, ", length =", length.sigma.adj2,'\n') W = 10000 ### CI for sigma using Weerahandi and s and s.ub r.sigma2 = rep(0,W) rub.sigma2 = rep(0,W) for (wb in 1:W) { wee = rnorm(1, 0, 1) r.sigma2[wb] = ( (wee+n/sqrt(2*n))*sqrt(2*n)/(s2*(n-1)) )^(-1) rub.sigma2[wb] = ( (wee+n/sqrt(2*n))*sqrt(2*n)/(s.ub^2*(n-1)) )^(-1) } l.score = max(0, quantile(r.sigma2, 0.025, na.rm=TRUE) ) u.score = quantile(r.sigma2, 0.975, na.rm=TRUE) length.score = sqrt(u.score)-sqrt(l.score) cat("CIgci.S =", sqrt(l.score), sqrt(u.score), ", length =", length.score,'\n') l.score.ub = max(0, quantile(rub.sigma2, 0.025, na.rm=TRUE) ) u.score.ub = quantile(rub.sigma2, 0.975, na.rm=TRUE) length.score.ub = sqrt(u.score.ub)-sqrt(l.score.ub) cat("CIgci.S.ub =", sqrt(l.score.ub), sqrt(u.score.ub ), ", length =", length.score.ub ,'\n') ### Estimated CV tau.ori = s/xbar tau.ub = s.ub/xbar var.tau.ub = (n*s2/(s2+n*xbar^2))*((n-1)/2)*(gamma((n-1)/2)/ gamma(n/2))^2 - (s/xbar)^2 tau.wp = sqrt((n-1)*s2/n )/xbar cat("tau.hat =", tau.ori*100, ", tau.hat.ub =", tau.ub*100, ", tau.hat.wp =", tau.wp*100,'\n') ### CI for mu l.mu = xbar - 1.96*s.ub/sqrt(n) u.mu = xbar + 1.96*s.ub/sqrt(n) ### CI for sigma using s.ub and wald method l.sigma = s.ub - 1.96*sqrt(var.s.ub) u.sigma = s.ub + 1.96*sqrt(var.s.ub) ### CI for sigma using original method in the book l.sigma.ori = sqrt(n-1)*s/sqrt(qchisq(1-alpha/2, df=n-1 )) u.sigma.ori = sqrt(n-1)*s/sqrt(qchisq(alpha/2, df=n-1 )) ### CI for sigma using adjusted original method (aa=1) l.sigma.adj = sqrt(n-1)*s.ub/sqrt(qchisq(1-alpha/2, df=n-1 )) u.sigma.adj = sqrt(n-1)*s.ub/sqrt(qchisq(alpha/2, df=n-1 )) l.sigma.adj2 = sqrt(n-aa)*s.ub/sqrt(qchisq(1-alpha/2, df=n-1 )) u.sigma.adj2 = sqrt(n-aa)*s.ub/sqrt(qchisq(alpha/2, df=n-1 )) ### CI_dz for tau from MOVER (Donner and zou) a = sqrt((n-1)/qchisq(1-alpha/2, df = n-1)) b = sqrt((n-1)/qchisq(alpha/2, df = n-1)) c = xbar^2 - 1.96^2*s2/n lower.tau.dz = (xbar - sqrt(max(0, xbar^2+a*c*(a-2))) )*s/c upper.tau.dz = (xbar + sqrt(max(0, xbar^2+b*c*(b-2))) )*s/c length.tau.dz = upper.tau.dz-lower.tau.dz cat("CIdz =", lower.tau.dz*100, upper.tau.dz*100, ", length =", length.tau.dz*100,'\n') #### CI_adj aa=2 for tau using adjusted original method 2 (with aa=2) lower.tau.adj2 = (s.ub*xbar - sqrt((s.ub*xbar)^2 - l.sigma.adj2*u.mu* (2*s.ub-l.sigma.adj2)*(2*xbar-u.mu)))/(u.mu*(2*xbar-u.mu)) upper.tau.adj2 = (s.ub*xbar + sqrt((s.ub*xbar)^2 - u.sigma.adj2*l.mu* (2*s.ub-u.sigma.adj2)*(2*xbar-l.mu)))/(l.mu*(2*xbar-l.mu)) length.tau.adj2 = upper.tau.adj2-lower.tau.adj2 cat("CImov.aa=2 =",lower.tau.adj2*100, upper.tau.adj2*100, ", length =", length.tau.adj2*100,'\n') ### CI for tau from Weerahandi using S and S.ub lower.tau.score = (s*xbar - sqrt((s*xbar)^2 - sqrt(l.score)*u.mu* (2*s-sqrt(l.score))*(2*xbar-u.mu)))/(u.mu*(2*xbar-u.mu)) upper.tau.score = (s*xbar + sqrt((s*xbar)^2 - sqrt(u.score)*l.mu* (2*s-sqrt(u.score))*(2*xbar-l.mu)))/(l.mu*(2*xbar-l.mu)) length.tau.score = upper.tau.score-lower.tau.score cat("CIgci.mov.gci.S =", lower.tau.score*100, upper.tau.score*100, ", length =", length.tau.score*100,'\n') lower.tau.score.ub = (s.ub*xbar - sqrt((s.ub*xbar)^2 - sqrt(l.score.ub)* u.mu*(2*s.ub-sqrt(l.score.ub))*(2*xbar-u.mu)))/(u.mu*(2*xbar-u.mu)) upper.tau.score.ub = (s.ub*xbar + sqrt((s.ub*xbar)^2 - sqrt(u.score.ub)* l.mu*(2*s.ub-sqrt(u.score.ub))*(2*xbar-l.mu)))/(l.mu*(2*xbar-l.mu)) length.tau.score.ub = upper.tau.score.ub-lower.tau.score.ub cat("CIgci.mov.gci.S.ub =",lower.tau.score.ub*100, upper.tau.score.ub*100, ", length =", length.tau.score.ub*100 ,'\n') #### CI from wararit k.tilda = sqrt((n-1)*s2/n )/xbar lower.tau.4 = k.tilda*( ( ((qchisq(1-alpha/2, df=n-1)+2)/(n))-1 )*k.tilda^2 + qchisq(1-alpha/2, df=n-1)/(n-1) )^(-0.5) upper.tau.4 = k.tilda*( ( ((qchisq(alpha/2, df=n-1)+2)/(n))-1 )*k.tilda^2 + qchisq(alpha/2, df=n-1)/(n-1) )^(-0.5) length.tau.4 = upper.tau.4-lower.tau.4 cat("CIwp =", lower.tau.4*100, upper.tau.4*100, ", length =", length.tau.4*100,'\n') } ########################### ### Put your data here ########################### x = c(1948407.81, 2010885.55 ,1411186.32, 1355634.83, 1332581.33, 1773235.72, 3368029.26, 3929424.56, 1584506.19, 1546560.88, 2388399.51, 866688.83, 496572.19, 1324523.97, 1749892.65, 904974.71, 1344246.07, 1659465.97, 822524.84, 1230819.33) ### Taking outliers out of the data extreme = quantile(x, prob = 0.75) + 3*IQR(x) out = boxplot(x)$out x = x[x <= extreme] n = length(x) ####################### ### Outputs are here! ####################### estimation(x, n)