c This program obtains a two-sided confidence c interval on a coefficient of variation for c data from a normal distribution. c c c Steve Verrill, 4/24/02 c Steve Verrill, 5/4/05 c c c implicit double precision (a-h,o-z) real*8 mean,sd,negb common /cvc/ estcv,sqrtn,df,ratio,alphad2,omad2 external flow,fup open(unit=7,file="dcdcvnorm.res",status="unknown") 5 write(6,10) write(7,10) 10 format(/,1x,'What is the sample mean?',/) read(5,*) mean write(7,*) mean if (mean .eq. 0.0d0) then write(6,12) write(7,12) 12 format(/,1x,'This program does not handle the', x ' case in which the sample mean equals 0.',/) go to 200 endif write(6,14) write(7,14) 14 format(/,1x,'What is the sample standard deviation?',/) read(5,*) sd write(7,*) sd if (sd .eq. 0.0d0) then write(6,16) write(7,16) 16 format(/,1x,'The confidence interval on the ', x 'coefficient of variation',/,1x,'is the single', x ' point 0.',/) go to 200 endif estcv = sd/mean write(6,18) estcv write(7,18) estcv 18 format(/,1x,'The estimated coefficient of variation', x' is ',e12.4,/) if (mean .lt. 0.0d0) then estcv = -estcv endif write(6,20) write(7,20) 20 format(/,1x,'What is n?',/) read(5,*) n write(7,*) n xn = n sqrtn = sqrt(xn) df = xn - 1.0d0 ratio = sqrtn/estcv 25 write(6,30) write(7,30) 30 format(/,1x,'What confidence level do you want?', x' (.95 for a 95% confidence level)',/) read(5,*) conflev write(7,*) conflev alphad2 = (1.0d0 - conflev)/2.0d0 omad2 = 1.0d0 - alphad2 if (n .gt. 3000) go to 1000 c Obtain a value for t_{n-1}^{-1}(1 - alphad2) call cdft(2,omad2,alphad2,tlim,df,ier,bound) if (ier .ne. 0) then write(6,32) write(7,32) 32 format(/,1x,'The ier value from a call', x ' to cdft was nonzero. Please contact',/, x 1x,'Steve Verrill at sverrill@fs.fed.us.',/) stop endif if (tlim .lt. 0.0d0 .or. tlim .gt. 100.0d0) then write(6,34) write(7,34) 34 format(/,1x,'It appears that', x ' t_{n-1}^{-1}(1 - alphad2) is less than 0',/, x 1x,'or greater than 100. As currently written,',/, x 1x,'this program cannot handle this sample size,',/, x 1x,'confidence level combination.',/) stop endif c c Obtain a lower confidence bound on cv (upper confidence bound c for a negative cv) c b1 = .001d0 c = 10.0d0 r = estcv re = .000001d0 ae = .00001d0 call dfzero(flow,b1,c,r,re,ae,iflag1) if (iflag1 .eq. 4 .and. b1 .lt. 5.0) then if (mean .gt. 0.0d0) then write(6,35) write(7,35) 35 format(/,1x,'The lower confidence bound on the COV', x /,1x,'lies between 0 and .001.',/) go to 100 else write(6,37) write(7,37) 37 format(/,1x,'The upper confidence bound on the COV', x /,1x,'lies between -.001 and 0.',/) go to 100 endif endif if (iflag1 .eq. 4 .and. b1 .gt. 5.0) then b2 = 10.0d0 c = 100.0d0 r = estcv re = .000001d0 ae = .001d0 call dfzero(flow,b2,c,r,re,ae,iflag2) if (iflag2 .eq. 4 .and. b2 .gt. 50.0d0) then if (mean .gt. 0.0d0) then write(6,45) write(7,45) 45 format(/,1x,'The lower confidence bound on the COV', x /,1x,'is greater than 100.',/) go to 200 else write(6,47) write(7,47) 47 format(/,1x,'The upper confidence bound on the COV', x /,1x,'is less than -100.',/) go to 200 endif else if (mean .gt. 0.0d0) then write(6,49) b2 write(7,49) b2 49 format(/,1x,'The lower confidence bound on the COV', x /,1x,'is ',e12.4,'.',/) go to 100 else negb = -b2 write(6,51) negb write(7,51) negb 51 format(/,1x,'The upper confidence bound on the COV', x /,1x,'is ',e12.4,'.',/) go to 100 endif endif endif c c iflag1 .ne. 4 c if (mean .gt. 0.0d0) then write(6,53) b1 write(7,53) b1 53 format(/,1x,'The lower confidence bound on the COV', x /,1x,'is ',e12.4,'.',/) go to 100 else negb = -b1 write(6,55) negb write(7,55) negb 55 format(/,1x,'The upper confidence bound on the COV', x /,1x,'is ',e12.4,'.',/) endif c c Obtain an upper confidence bound on cv (lower confidence bound c for a negative cv) c c c If the tlim obtained above is greater than c or equal to sqrt(n)/estcv, then the upper confidence c bound is positive infinity. c 100 if (tlim .ge. ratio) then if (mean .gt. 0.0d0) then write(6,110) write(7,110) 110 format(/,1x,'The upper confidence bound on the COV', x ' is positive infinity.',/) go to 200 else write(6,112) write(7,112) 112 format(/,1x,'The lower confidence bound on the COV', x ' is negative infinity.',/) go to 200 endif endif b1 = .001d0 c = 10.0d0 r = estcv re = .000001d0 ae = .00001d0 call dfzero(fup,b1,c,r,re,ae,iflag1) if (iflag1 .eq. 4 .and. b1 .lt. 5.0) then if (mean .gt. 0.0d0) then write(6,135) write(7,135) 135 format(/,1x,'The upper confidence bound on the COV', x /,1x,'lies between 0 and .001.',/) go to 200 else write(6,137) write(7,137) 137 format(/,1x,'The lower confidence bound on the COV', x /,1x,'lies between -.001 and 0.',/) go to 200 endif endif if (iflag1 .eq. 4 .and. b1 .gt. 5.0) then b2 = 10.0d0 c = 100.0d0 r = estcv re = .000001d0 ae = .001d0 call dfzero(fup,b2,c,r,re,ae,iflag2) if (iflag2 .eq. 4 .and. b2 .gt. 50.0d0) then if (mean .gt. 0.0d0) then write(6,145) write(7,145) 145 format(/,1x,'The upper confidence bound on the COV', x /,1x,'is greater than 100.',/) go to 200 else write(6,147) write(7,147) 147 format(/,1x,'The lower confidence bound on the COV', x /,1x,'is less than -100.',/) go to 200 endif else if (mean .gt. 0.0d0) then write(6,149) b2 write(7,149) b2 149 format(/,1x,'The upper confidence bound on the COV', x /,1x,'is ',e12.4,'.',/) go to 200 else negb = -b2 write(6,151) negb write(7,151) negb 151 format(/,1x,'The lower confidence bound on the COV', x /,1x,'is ',e12.4,'.',/) go to 200 endif endif endif c c iflag1 .ne. 4 c if (mean .gt. 0.0d0) then write(6,153) b1 write(7,153) b1 153 format(/,1x,'The upper confidence bound on the COV', x /,1x,'is ',e12.4,'.',/) go to 200 else negb = -b1 write(6,155) negb write(7,155) negb 155 format(/,1x,'The lower confidence bound on the COV', x /,1x,'is ',e12.4,'.',/) endif go to 200 c Use the McKay-Vangel approach for n greater tha 3000 1000 xk = estcv xk2 = xk*xk c c Lambda_4 c call cdfchi(2,omad2,alphad2,u1,df,istatus,bound) if (istatus .ne. 0) then write(6,1010) write(7,1010) 1010 format(/,1x,'The ier value from a call', x ' to cdfchi was nonzero. Please contact',/, x 1x,'Steve Verrill at sverrill@fs.fed.us.',/) stop endif call cdfchi(2,alphad2,omad2,u2,df,istatus,bound) if (istatus .ne. 0) then write(6,1010) stop endif arg = ((u1+2)/(df + 1.0d0) - 1)*xk2 + u1/df blm = xk/sqrt(arg) if (mean .gt. 0.0d0) then write(6,1020) blm write(7,1020) blm 1020 format(/,1x,'The lower confidence bound on the COV', x /,1x,'is ',e12.4,'.',/) go to 1100 else blmm = -blm write(6,1030) blmm write(7,1030) blmm 1030 format(/,1x,'The upper confidence bound on the COV', x /,1x,'is ',e12.4,'.',/) endif 1100 arg = ((u2+2)/(df + 1.0d0) - 1)*xk2 + u2/df if (arg .gt. 0.0d0) then bum = xk/sqrt(arg) if (mean .gt. 0.0d0) then write(6,1120) bum write(7,1120) bum 1120 format(/,1x,'The upper confidence bound on the COV', x /,1x,'is ',e12.4,'.',/) go to 200 else bumm = -bum write(6,1130) bumm write(7,1130) bumm 1130 format(/,1x,'The lower confidence bound on the COV', x /,1x,'is ',e12.4,'.',/) go to 200 endif else if (mean .gt. 0.0d0) then write(6,1220) write(7,1220) 1220 format(/,1x,'The upper confidence bound on the COV', x /,1x,'is plus infinity.',/) go to 200 else write(6,1230) write(7,1230) 1230 format(/,1x,'The lower confidence bound on the COV', x /,1x,'is negative infinity.',/) go to 200 endif endif 200 write(6,210) 210 format(/,1x,'Another calculation? 0 - no 1 - yes',/) read(5,*) iflag if (iflag .eq. 1) go to 5 stop end double precision function fup(cv) implicit double precision (a-h,o-z) common /cvc/ estcv,sqrtn,df,ratio,alphad2,omad2 xncp = sqrtn/cv call cdftnc(1,p,q,ratio,df,xncp,ier,bound) if (ier .ne. 0) then write(6,10) write(7,10) 10 format(/,1x,'The ier value from a call', x ' to cdftnc was nonzero. Please contact',/, x 1x,'Steve Verrill at sverrill@fs.fed.us.',/) stop endif fup = omad2 - p return end double precision function flow(cv) implicit double precision (a-h,o-z) common /cvc/ estcv,sqrtn,df,ratio,alphad2,omad2 xncp = sqrtn/cv call cdftnc(1,p,q,ratio,df,xncp,ier,bound) if (ier .ne. 0) then write(6,10) write(7,10) 10 format(/,1x,'The ier value from a call', x ' to cdftnc was nonzero. Please contact',/, x 1x,'Steve Verrill at sverrill@fs.fed.us.',/) stop endif flow = alphad2 - p return end