*! 1.1.0 NJC 25 November 2016
*! 1.0.0 NJC 2 February 2009
program diptest, sort rclass
	version 9
	syntax varname(numeric) [if] [in] [, by(varlist) Reps(int 10000) * ]

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

		tempvar group guse runif
		tempname dip P low high mean sim

		if "`by'" == "" {
			tempvar by
			gen byte `by' = 1
		}

		egen `group' = group(`by') if `touse'
		su `group' if `touse', meanonly
		local ng = r(max)

		gen byte `guse' = .
		gen double `runif' = .
		sort `touse' `group' `varlist'

		tokenize n dip P low high mean
		forval j = 1/6 {
			tempvar show`j'
  			gen `show`j'' = .
			char `show`j''[varname] "``j''"
			local showlist `showlist' `show`j''
		}

		format %6.4f `show2' `show3'

		forval g = 1/`ng' {
			replace `guse' = `touse' & (`group' == `g')
	        count if `guse'
			return scalar n_`g' = r(N)
			replace `show1' = r(N) if `group' == `g'

			mata: ///
		diptest("`varlist'", "`guse'", "`dip'", "`low'", "`high'", "`mean'" )
			return scalar dip_`g' = scalar(`dip')
			return scalar low_`g' = scalar(`low')
			return scalar high_`g' = scalar(`high')
			return scalar mean_`g' = scalar(`mean')
			replace `show2' = scalar(`dip') if `group' == `g'
			replace `show4' = scalar(`low') if `group' == `g'
			replace `show5' = scalar(`high') if `group' == `g'
			replace `show6' = scalar(`mean') if `group' == `g'

			local freq = 0
			forval i = 1/`reps' {
				replace `runif' = uniform()
				mata: ///
		diptest("`runif'", "`guse'", "`sim'", "`low'", "`high'", "`mean'")
				local freq = `freq' + (scalar(`sim') > scalar(`dip'))
			}

			scalar `P' = `freq'/`reps'
			return scalar P_`g' = scalar(`P')
			replace `show3' = scalar(`P') if `group' == `g'
		} // loop over groups

		bysort `touse' `group': replace `touse' = `touse' & (_n == _N)
	}

	if `ng' == 1 list `showlist' if `touse', subvarname noobs `options'
	else list `by' `showlist' if `touse', subvarname noobs `options'

end

// diptest NJC 2 February 2009
// based on Fortran code in
// Hartigan, P.M. 1985. Algorithm AS 217: Computation of the dip
// statistic to test for unimodality. Applied Statistics 34: 320-325.
// and C code in
// Maechler, Martin. 2003. diptest 0.25-1. http://www.r-project.org/

mata:

void diptest(string scalar varname,
	     string scalar tousename,
	     string scalar dipname,
	     string scalar lowname,
	     string scalar highname,
	     string scalar meanname
	) {

	real colvector x, mn, mj, lcm, gcm
	real scalar n, high, low, gcmi, gcmi1, gcmix,
		ic, icx, icv, ig, ih, iv, ix, j, jb, je, jr,
		k, kb, ke, kr, lcm1, lcmiv, lcmiv1, mjk, mjmjk, mnj, mnmnj
	real scalar dx, dip, d, dl, du, dipnew
	real scalar temp, C, t

	x = st_data(., varname, tousename)

	if (rows(x) == 0 | cols(x) == 0) {
		errprintf("no data")
		return
	}

	_sort(x, 1)
	n = rows(x)

	if (n <= 3 | x[1] == x[n]) {
		st_numscalar(dipname, 0)
		st_numscalar(lowname, x[1])
		st_numscalar(highname, x[n])
		st_numscalar(meanname, mean(x))
		return
	}

	low = 1 ; high = n ; dip = 1 / n

	mn = 1 \ (1 :: n - 1)
	for(j = 2; j <= n; j++) {
		while (1) {
	        	mnj = mn[j]
			mnmnj = mn[mnj]
			if (mnj == 1 |
(mnj - mnmnj) * (x[j] - x[mnj]) < (j - mnj) * (x[mnj] - x[mnmnj])) break
			mn[j] = mnmnj
		}
	}

	mj = mn; mj[n] = n
	for(j = 1; j < n; j++) {
		k = n - j
		mj[k] = k + 1
		while (1) {
                	mjk = mj[k]
			mjmjk = mj[mjk]
			if (mjk == n |
(mjk - mjmjk) * (x[k] - x[mjk]) < (k - mjk) * (x[mjk] - x[mjmjk])) break
			mj[k] = mjmjk
		}
	}

    p1: gcm = mn; gcm[1] = high
	for(ic = 1; gcm[ic] > low; ic++) {
		gcm[ic + 1] = mn[gcm[ic]]
	}
	icx = ic

	lcm = mn; lcm[1] = low
	for(ic = 1; lcm[ic] < high; ic++) {
		lcm[ic + 1] = mj[lcm[ic]]
	}
	icv = ic

	ig = icx; ih = icv ; ix = icx - 1; iv = 2

	d = 0
	if (icx != 2 | icv != 2) goto p2
	d = 1/n
	goto p3

    p2: gcmix = gcm[ix]
	lcmiv = lcm[iv]
	if (gcmix <= lcmiv) {
		lcmiv1 = lcm[iv - 1]
			dx = (x[gcmix] - x[lcmiv1]) * (lcmiv - lcmiv1) /
            (n * (x[lcmiv] - x[lcmiv1])) - (gcmix - lcmiv1 - 1) / n

			// bug fix 25 Nov 2016: was
		    // dx = (x[gcmix] - x[lcmiv1] * (lcmiv - lcmiv1)) /
			// (n * (x[lcmiv] - x[lcmiv1])) - (gcmix - lcmiv1 - 1) / n
		ix--
		if (dx >= d) {
			d = dx
			ig = ix + 1
			ih = iv
		}
	}
	else {
		lcmiv = lcm[iv]
		gcmi = gcm[ix]
		gcmi1 = gcm[ix + 1]
		dx = (lcmiv - gcmi1 + 1) / n -
			(x[lcmiv] - x[gcmi1]) * (gcmi - gcmi1) / (n * (x[gcmi] - x[gcmi1]))
		iv++
		if (dx >= d) {
			d = dx
			ig = ix + 1
			ih = iv - 1
		}
	}

	if (ix < 1) ix = 1
	if (iv > icv) iv = icv
	if (gcm[ix] != lcm[iv]) goto p2

    p3:	if (d < dip) goto done

	dl = 0
	if (ig != icx) {
		for(j = ig; j < icx; j++) {
			temp = 1/n
			jb = gcm[j + 1]
			je = gcm[j]
			if ((je - jb > 1) & (x[je] != x[jb])) {
				C = (je - jb) / (n * (x[je] - x[jb]))

				for(jr = jb; jr <= je; jr++) {
					t = (jr - jb + 1) / n - (x[jr] - x[jb]) * C
					if (t > temp) temp = t
				}
			}
	                if (dl < temp) dl = temp
		}
	}

        du = 0
        if (ih != icv) {
		for(k = ih; k < icv; k++) {
			temp = 1/n
			kb = lcm[k]
			ke = lcm[k + 1]
			if ((ke - kb > 1) & (x[ke] != x[kb])) {
				C = (ke - kb) / (n * (x[ke] - x[kb]))
				for (kr = kb; kr <= ke; kr++) {
					t = (x[kr] - x[kb]) * C - (kr - kb - 1) / n
					if (t > temp) temp = t
				}
			}
            	        if (du < temp) du = temp
		}
	}

        dipnew = du > dl ? du : dl
	if (dip < dipnew) dip = dipnew

	if (low == gcm[ig] & high == lcm[ih]) {
		// do nothing
	}
	else {
		low  = gcm[ig]; high = lcm[ih]
		goto p1
	}

  done: st_numscalar(dipname, dip / 2)
        st_numscalar(lowname, x[low])
	st_numscalar(highname, x[high])
	st_numscalar(meanname, mean(x[low::high]))

}

end