*! NJC 1.0.0 16 Nov 2006
* Pareto dot plot
* Wilkinson, Leland. 2006. Revising the Pareto chart.
* American Statistician 60(4): 332-334.
program pdplot, sort
	version 9
	syntax varname [if] [in] [fweight/] ///
	[, DOTSonly nreps(int 10000) Horizontal AIopts(str asis) ///
	Level(int `c(level)') * ]

	quietly {
		marksample touse, strok
		count if `touse'
		if r(N) == 0 error 2000

		tempvar freq tag rank order
		tempname ranklbl low high

		if "`exp'" == "" local exp = 1

		// frequencies
		bysort `touse' `varlist' : gen `freq' = sum(`exp')
		by `touse' `varlist' : replace `freq' = `freq'[_N]

		// categories
		bysort `touse' `freq' `varlist' : ///
			gen byte `tag' = _n == 1 & `touse'
		count if `tag'
		local ncats = r(N)

		su `freq' if `tag', meanonly
		local nvals = r(sum)

		// ranks
		gen `rank' = `ncats' - sum(`tag') + 1
		local label : var label `varlist'
		if `"`label'"' == "" local label "`varlist'"
		label var `rank' `"`label'"'

		// label ranks
		gen long `order' = _n
		capture confirm string var `varlist'
		local isnum = _rc

		forval i = 1/`ncats' {
			su `order' if `rank' == `i', meanonly
			local value = `varlist'[r(min)]
			if `isnum' local value : label (`varlist') `value'
			label def `ranklbl' `i' `"`value'"', modify
		}

		if "`dotsonly'" == "" {
			// 95% acceptance intervals
			mata: ///
		work(`nreps',`ncats',`nvals',`level',"`c(seed)'","`low'","`high'")

			tempvar lo hi
			gen `lo' = .
			gen `hi' = .
			char `lo'[varname]   "Lower"
			char `hi'[varname]   "Upper"

			forval i = 1/`ncats' {
				local j = `ncats' - `i' + 1
				replace `hi' = `high'[`j',1] if `rank' == `i'
				replace `lo' = `low'[`j',1] if `rank' == `i'
			}
		}

		sort `touse' `tag' `rank'

		// percent scale too
		local pcmax = 100 * `freq'[_N - `ncats' + 1] / `nvals'
		local step = cond(`pcmax' < 25, 5, 10)
		local pcmax = `step' * ceil(`pcmax'/`step')

		forval v = 0(`step')`pcmax' {
			local val = `v' * `nvals'/100
			local percents `"`percents' `val' "`v'""'
		}
	}

	// list
	char `freq'[varname] "Freq."
	char `rank'[varname] "Rank"
	list `varlist' `rank' `freq' `lo' `hi' if `tag', ///
		subvarname noobs sep(0)

	if "`dotsonly'" == "" {
		di as res "  (`level'% acceptance intervals from `nreps' random samples)"
	}

	// graph
	label val `rank' `ranklbl'

	if "`dotsonly'" == "" {
		#delimit ;
		local bars  rbar `lo' `hi' `rank' if `tag',
    		`horizontal'
    		bcolor(none)
    		barw(0.2)
                note(`level'% acceptance intervals)
		`aiopts'
		|| ;
		#delimit cr
	}

	if "`horizontal'" != "" {
		#delimit ;
		twoway `bars'
		scatter `rank' `freq' if `tag',
		xla(, ang(h))
		xaxis(1 2)
		xla(`percents', ang(h) axis(1))
		xtitle(Percent, axis(1))
		xtitle(Frequency, axis(2))
		yaxis(1 2)
		yla(1/`ncats', ang(h) axis(1))
		yla(1/`ncats', ang(h) valuelabel noticks axis(2))
		ytitle(Rank, axis(1)) ysc(reverse)
		legend(off)
		`options' ;
		#delimit cr
	}
	else {
		#delimit ;
		twoway `bars'
		scatter `freq' `rank' if `tag',
		yaxis(1 2)
		yla(`percents', ang(h) axis(1))
		yla(, ang(h) axis(2))
		ytitle(Percent, axis(1))
		ytitle(Frequency, axis(2))
		xaxis(1 2)
		xla(1/`ncats', valuelabel noticks axis(2))
		xla(1/`ncats', axis(1))
		xtitle(Rank, axis(1)) ms(Oh)
		legend(off)
		`options' ;
		#delimit cr
	}
end

mata:

void work(real scalar nreps,
          real scalar ncats,
	  real scalar nvals,
	  real scalar level,
	  string scalar seed,
	  string scalar lowname,
	  string scalar highname)
{
	real colvector sample, low, high, counts
	real matrix allcounts
	real scalar i, j, nlow, nhigh

	counts = J(ncats, 1, 0)
	allcounts = J(nreps, ncats, .)
	low = high = J(ncats, 1, 0)
	uniformseed(seed)

	for(i = 1; i <= nreps; i++) {
		sample = ceil(ncats * uniform(nvals, 1))
		for(j = 1; j <= ncats; j++) {
			counts[j,1] = sum(sample :== j)
		}
		counts = sort(counts, 1)
		allcounts[i,] = counts'
	}

	// explanation below
	level = (100 - level) / 200
	nlow = floor(0.5 + nreps * level)
	nhigh = ceil(0.5 + nreps * (1 - level))

	for(j = 1; j <= ncats; j++) {
		allcounts = sort(allcounts, j)
		low[j,1] = allcounts[nlow, j]
		high[j,1] = allcounts[nhigh, j]
	}

	st_matrix(lowname, low)
	st_matrix(highname, high)
}

end

/*

Given level as a percent: say 95 for example.

We work with alpha / 2 = (100 - level) / 200: in this example 0.025.

We seek integer solutions i to
	alpha / 2     = (i - 0.5) / n
	1 - alpha / 2 = (i - 0.5) / n
which will be integers close to
	0.5 + n * alpha / 2
	0.5 + n * (1 - alpha / 2)
Arbitrarily, we round down the lower value using floor() and round up
the upper value using ceil().

*/