* Stata 3.0 version 1.0.0 JPR 04.03.92.
* CENTILES - nonparametric or normal-based centile estimation
* with confidence limits. Weights not yet implemented.
program define centile
	version 3.0
	local varlist "opt ex min(1)"
	local if "opt"
	local in "opt"
	local options /*
		*/ "CCi Centile(string) Level(integer $S_level) Meansd Normal"
	parse "`*'"
	if "`centile'"=="" { local centile 50 }
/*
	Parse `centile' and count the requested centiles, placing them
	into c1, c2, ....
*/
	local nc 0
	if "`centile'" ~= "" {
		parse "`centile'", parse(" ,")
		while "`1'" ~= "" {
			if "`1'" == "," { 
				mac shift 
			}
			confirm number `1'
			local nc = `nc'+1
			local c`nc' `1'
			local cj "c`nc'"
			mac shift
		}
	}
	local tl1 "    Obs "
	local ttl "Percentile"
	if "`meansd'"=="" { 
		if "`normal'"=="" { 
			if "`cci'"~="" {
				di in gr _n _col(52) "-- Binomial Exact --"
			}
			else {
				di in gr _n _col(52) "-- Binom. Interp. --"
			}
		}
		else {
			di in gr _n _col(32) "-- Normal, based on observed centiles --"
		}
	}
	else {
		di in gr _n _col(32) "-- Normal, based on mean and std. dev.--"
	}
	#delimit ;
	di in gr 
"Variable | `tl1' `ttl'      Centile        [`level'% Conf. Interval]" 
	_n "---------+" _dup(61) "-" ;
	#delimit cr

	local anymark 0
	local alpha2 = (100-`level')/200
	local zalpha2 = -invnorm(`alpha2')
	tempvar YVAR
	qui gen `YVAR' = .
/*
	Run through varlist
*/
	parse "`varlist'", parse(" ")
	local vl
	while "`1'" != "" {
		capt conf str v `1'
		if _rc { 
			local vl "`vl' `1'"
		}
		mac shift 
	}

	local nobs 0			/* in case loop aborted */
	parse "`vl'", parse(" ")
	while "`1'" ~= "" {
		local yvar "`1'"
		qui replace `YVAR' = `yvar' `if' `in'
		sort `YVAR'
		qui sum `YVAR'
		local nobs = _result(1)
		local mean = _result(3)
		local sd = sqrt(_result(4))
/*
	Formatting
*/
		local skip = 8 - length("`yvar'")
		local fmt : format `yvar'
		if substr("`fmt'",-1,1)=="f" { 
			local ofmt="%9."+substr("`fmt'",-2,2)
		}
		else	local ofmt "%9.0g"
/*
	Calc required centile(s)
*/
		local j 1
		local s 7
		while `j' <= `nc' {
			local mark ""
			local cj "c`j'"
			local quant = ``cj''/100
			if "`meansd'" ~= "" {
/*
	Normal distribution (parametric estimates)
*/
				local z = invnorm(`quant')
				local centil = `mean'+`z'*`sd'
				local se = `sd'*sqrt(1/`nobs'+`z'^2/(2*(`nobs'-1)))
				local cLOWER = `centil'-`zalpha2'*`se'
				local cUPPER = `centil'+`zalpha2'*`se'
			}
			else {
/*
	If `normal' is not set, calc centile and exact
	nonparametric confidence limits (for example,
	see Gardner and Altman 1989, pp 72-74),
	interpolating when not at ends of distribution.
	An iterative process is required for each limit.
	As starting values, find ranks for lower and
	upper limits using a normal approximation.

	If `normal' is set, use normal distribution formula
	for variance, (10.29) in Kendall & Stuart (1969) p 237.

	`alpha2' is for example .025 for a 95% CI.
*/
				local frac1 = (`nobs'+1)*`quant'
* di "Central fraction and rank = " `quant',`frac1'
				local i1 = int(`frac1')
				local frac1 = `frac1'-`i1'
				if `i1' >= `nobs' {
					local centil = `YVAR'[`nobs']
				}
				else if `i1' < 1 {
					local centil = `YVAR'[1]
				}
				else {
					local centil = `YVAR'[`i1']+`frac1'*(`YVAR'[`i1'+1]-`YVAR'[`i1'])
				}
				if "`normal'" == "" {
					local nq = `nobs'*`quant'
					local z = sqrt(`nq'*(1-`quant'))*`zalpha2'
					local rzLOW = int(.5+`nq'-`z')
					local rzHIGH = 1+int(.5+`nq'+`z')
* di "lower and upper approx ranks = " `rzLOW',`rzHIGH'
					local r1 `rzHIGH'
					if `r1' > `nobs'+1 { local r1 = `nobs'+1 }
					local r0 = `r1'-1
					local p0 = Binomial(`nobs',`r0',`quant')
					local p1 = Binomial(`nobs',`r1',`quant')
					local done 0
					while ~`done' {
						if `p0' > `alpha2' {
							if `p1' <= `alpha2' {
								local done 1
							}
							else {
								local r0 = `r1'
								local p0 = `p1'
								local r1 = `r1'+1
								local p1 = Binomial(`nobs',`r1',`quant')
							}
						}
						else if `p0' == `alpha2' {
							local r1 = `r0'
							local p1 = `p0'
							local done 1
						}
						else {
							local r1 = `r0'
							local p1 = `p0'
							local r0 = `r0'-1
							local p0 = Binomial(`nobs',`r0',`quant')
						}
					}
* di "Upper r0, p0, r1, p1, interp =" `r0',`p0',`r1',`p1',(`p0'-`alpha2')/(`p0'-`p1')
/*
	Upper limit. Interpolate between r0 and r1, r1 being
	conservative. Note that p0>=p1 (upper tail areas).
*/
					if `r0' >= `nobs' {
						local cUPPER = `YVAR'[`nobs']
						local mark "*"
						local anymark 1
					}
					else if `r0' < 1 {
						local cUPPER = `YVAR'[1]
						local mark "*"
						local anymark 1
					}
					else {
						if "`cci'" == "" {
							local cUPPER = `YVAR'[`r0']		/*
							*/	+((`p0'-`alpha2')/(`p0'-`p1'))	/*
							*/	*(`YVAR'[`r1']-`YVAR'[`r0'])
						}
						else {
							local cUPPER = `YVAR'[`r1']
						}
					}
					local r1 `rzLOW'
					if `r1' < 0 { local r1 0 }
					local r0 = `r1'-1
					local p0 = 1-Binomial(`nobs',`r0'+1,`quant')
					local p1 = 1-Binomial(`nobs',`r1'+1,`quant')
					local done 0
					while ~`done' {
						if `p1' > `alpha2' {
							if `p0' <= `alpha2' {
								local done 1
							}
							else {
								local r1 = `r0'
								local p1 = `p0'
								local r0 = `r0'-1
								local p0 = 1-Binomial(`nobs',`r0'+1,`quant')
							}
						}
						else if `p0' == `alpha2' {
							local r0 = `r1'
							local p0 = `p1'
							local done 1
						}
						else {
							local r0 = `r1'
							local p0 = `p1'
							local r1 = `r1'+1
							local p1 = 1-Binomial(`nobs',`r1'+1,`quant')
						}
					}
/*
	Interpolate between r1 and r1+1, r1 being conservative.
	Note that p0<=p1 (both are lower tail areas).
*/
					if `r1' >= `nobs' {
						local cLOWER = `YVAR'[`nobs']
						local mark "*"
						local anymark 1
					}
					else if `r1' < 1 {
						local cLOWER = `YVAR'[1]
						local mark "*"
						local anymark 1
					}
					else {
						if "`cci'" == "" {
							local cLOWER = `YVAR'[`r1'] /*
							*/ +((`alpha2'-`p0')/(`p1'-`p0')) /*
							*/ *(`YVAR'[`r1'+1]-`YVAR'[`r1'])
						}
						else {
							local cLOWER = `YVAR'[`r1']
						}
					}
* di "Lower r0, p0, r1, p1, interp =" `r0',`p0',`r1',`p1',(`alpha2'-`p0')/(`p1'-`p0')
* di "Var = `yvar', centile `cj'=" `centil' ", CI = " `cLOWER',`cUPPER'
				}
				else {
/*
	Normal distribution, observed centiles
*/
					local dens = exp(-0.5*((`centil'-`mean')/`sd')^2)/(`sd'*sqrt(2*_pi))
					local se = sqrt(`quant'*(1-`quant')/`nobs')/`dens'
					local cLOWER = `centil'-`zalpha2'*`se'
					local cUPPER = `centil'+`zalpha2'*`se'
				}
			}
			if `j' == 1 {
				di in gr _skip(`skip') "`yvar' |" _col(10) /*
				*/ in yel %8.0f `nobs'		/*
		 		*/ _col(23) %7.0g ``cj''	/*
		 		*/ _col(35) `ofmt' `centil'	/*
		 		*/ _col(51) `ofmt' `cLOWER'	/*
		 		*/ _col(63) `ofmt' `cUPPER' in gr "`mark'"
			}
			else {
				di in gr "         |"		/*
		 		*/ in yel _col(23) %7.0g ``cj''	/*
		 		*/ _col(35) `ofmt' `centil'	/*
		 		*/ _col(51) `ofmt' `cLOWER'	/*
		 		*/ _col(63) `ofmt' `cUPPER' in gr "`mark'"
			}
/*
	Store centile in S_ macros, starting at 7.
*/
			mac def S_`s' `centil'
			local j = `j'+1
			local s = `s'+1
		}
		mac shift
	}
	if "`anymark'" == "1" {
		di in gr _n /*
*/ "`mark' Lower (upper) confidence limit held at minimum (maximum) of sample"
	}
/*
	Store quantities at final point; S_4 is duplicate of S_[`nc'+6].
*/
	mac def S_1 `nobs'
	mac def S_2 `nc'
	mac def S_3 ``cj''
	mac def S_4 `centil'
	mac def S_5 `cLOWER'
	mac def S_6 `cUPPER'
end