*!IHSU 1.0.3 03 January 2013
* IHSU 1.0.2 26 December 2012
* IHSU 1.0.1 17 November 2012
* NJC 2.0.1 2 April 2004
* NJC 2.0.0 21 January 2004
* NJC 1.1.1 16 December 1998
* NJC 1.1.0 26 October 1996
* nonparametric density estimation for circular data: Fisher 1993 pp.24-27
* Modified initially 17 November 2012 by Isaias Hazarmabeth Salgado-Ugarte
* This program calculates kernel density estimator of a series of values
* on a circular scale according to the weight functions coded as follows:
* 1 = Uniform
* 2 = Triangle
* 3 = Epanechnikov
* 4 = Quartic (Biweight)
* 5 = Triweight
* 6 = Gaussian
* 7 = Cosinus
* This version counts and estimates modes and antimodes
* and draws a linear or a circular graph
* Updated by IHSU & VMSQ 07 March 2021
program circkden
version 11.0
syntax varname(numeric) [if] [in] ///
[, H(real 30) Kc(integer 4) NPoints(integer 0) NUMOdes MOdes NUAMOdes AMOdes ///
NOGraph CIRCGph Rval(real 1) Fr(real 1) GS(real 1) GEN(str) PLOT(str asis) * ]
marksample touse
qui count if `touse'
if r(N) == 0 error 2000
local nd=r(N)
preserve
if "`gen'" != "" {
parse "`gen'", parse(" ")
confirm new variable `1'
confirm new variable `2'
}
tempvar deg d e f g cosdeg sindeg denp cosdenp sindenp
qui {
if `npoints'==0 {
if _N < 360 {
gen `deg' = (_n -1) / (_N-1) * 360
}
else gen `deg' = _n in 1/360
*noi li `deg'
}
else {
if `npoints'< r(N) {
di as error "The number of points must be >= " r(N)
exit
}
set obs `npoints'
count if `touse'
*noi di `npoints'
gen `deg' = (_n -1) / (`npoints'-1) * 360
}
*noi li `deg'
if `"`:variable label `varlist''"' != "" {
label var `deg' `"`:variable label `varlist''"'
}
else label var `deg' "`varlist'"
*local n = min(_N, 360)
local n = `nd'
gen `d' = .
gen `e' = .
gen `f' = .
if `kc' == 1 {
local klab = "(Uniform)"
}
else if `kc' == 2 {
local klab = "(Triangular)"
}
else if `kc' == 3 {
local klab = "(Epanechnikov)"
}
else if `kc' == 4 {
local klab = "(Biweight)"
}
else if `kc' == 5 {
local klab = "(Triweight)"
}
else if `kc' == 6 {
local klab = "(Gaussian)"
}
else if `kc' == 7 {
local klab = "(Cosine)"
}
local hlab=string(round(`h',.0001),"%8.2f")
label var `f' ///
"Density"
if "`modes'"~="" & "`numodes'"=="" {
di in red "you must include the 'numodes' option"
exit
}
if "`amodes'"~="" & "`nuamodes'"=="" {
di in red "you must include the 'nuamodes' option"
exit
}
local i = 1
while `i' <= _N {
replace `d' = abs(`deg'[`i'] - `varlist') if `touse'
replace `e' = min(`d' , 360 - `d')
if `kc'==1 {
replace `e'= cond(`e' >= `h', 0, 0.5)
}
else if `kc'==2 {
replace `e'= cond(`e' >= `h', 0, (1-abs(`e'/`h')))
}
else if `kc'==3 {
replace `e'= cond(`e' >=`h', 0, ((3/4)*(1-`e'^2/`h'^2)))
}
else if `kc'==4 {
replace `e'= cond(`e' >=`h',0, (15/16)*((1-`e'^2/`h'^2))^2)
}
else if `kc'==5 {
replace `e'=cond(`e' >=`h',0,((35/32)*(1-`e'^2/`h'^2)^3))
}
else if `kc'==6 {
replace `e'=(1/(sqrt(2*_pi)))*exp(-0.5*`e'^2/`h'^2)
}
else if `kc'==7{
replace `e'=cond(`e' >=`h',0,(_pi/4)*cos((_pi/2)*`e'/`h'))
}
else {
di as error "The kernel code must be an integer from 1 to 7"
exit
}
su `e', meanonly
replace `f' = r(sum) / (`n' * `h') in `i'
local i = `i' + 1
replace `f' = . if `deg' == .
}
if "`nograph'"=="" {
if "`circgph'"=="" {
if `"`subtitle'"' == "" {
local subtitle "sub("Circular kernel `klab' density estimate, ///
{it:h} = `h'{°ree}", size(medium))"
}
twoway scatter `f' `deg', ///
xla(0 "0" 90 "90" 180 "180" 270 "270" 360 "360") yla(, ang(h)) c(l) ///
`subtitle' `options' ///
|| `plot'
}
else {
gen `cosdeg' = cos(`deg'*_pi/180)
gen `sindeg' = sin(`deg'*_pi/180)
sum `f'
gen `g' = `f'/r(max)
gen `denp' = `rval'*(1 + _pi*`g')^.5 - `rval'
gen `cosdenp' = `cosdeg'*(1 + `denp'*`fr')
gen `sindenp' = `sindeg'*(1 + `denp'*`fr')
if `"`subtitle'"' == "" {
local subtitle "sub("Circular kernel `klab' density estimate, ///
{it:h} = `h'{°ree}", pos(6) size(medium))"
}
else local subtitle `"sub(`subtitle')"'
local size = 1.1 + `gs'
scatter `cosdeg' `sindeg', ms(i) c(l) || ///
scatter `cosdenp' `sindenp', ms(i) aspect(1) c(l) legend(off) yline(0) ///
xline(0) ysc(r(-`size' `size') off fill) xsc(r(-`size' `size') off fill) ylab(, nogrid) plotregion(margin(zero) style(none)) ///
`subtitle' `options' ///
|| `plot'
}
}
if "`numodes'"~="" {
tempvar difvar inmo sumo
gen `difvar'=`f'[_n+1] - `f'[_n] if `deg' != .
gen `inmo' = 0 if `deg' != .
replace `inmo'=1 if `difvar'[_n]>=0 & `difvar'[_n+1] < 0
gen `sumo' = sum(`inmo')
local numo= `sumo'[_N]
noi di as txt _newline " Number of modes = " as res `numo'
}
if "`modes'"~="" {
tempvar modes
gen `modes'=.
replace `modes'=`deg' if `inmo'[_n-1]==1
sort `modes'
local i = 1
noi di as txt _newline _dup(75) "_"
local title " Modes in circular KDE"
count if `deg' !=.
noi di as txt "`title', bw = " as res `h' as txt ", Ker = " as res `kc' as txt ", deg points = " as res r(N)
*noi di as txt "`title', bw = `h', Ker = `kc', npoints = `n'"
noi di as txt _dup(75) "-"
while `i'<`numo'+1 {
noi di as txt " Mode ( " as res %4.0f `i' as txt " ) = " as res %12.4f `modes'[`i']
local i = `i'+1
}
noi di as txt _dup(75) "_"
sort `deg'
}
if "`nuamodes'"~="" {
tempvar difvar inamo suamo
gen `difvar'=`f'[_n+1] - `f'[_n] if `deg' != .
gen `inamo' = 0 if `deg' != .
replace `inamo'=1 if `difvar'[_n]<=0 & `difvar'[_n+1] > 0
gen `suamo' = sum(`inamo')
local nuamo= `suamo'[_N]
noi di as txt _newline " Number of antimodes = " as res `nuamo'
}
if "`amodes'"~="" {
tempvar amodes
gen `amodes'=.
replace `amodes'=`deg' if `inamo'[_n-1]==1
sort `amodes'
local i = 1
noi di as txt _newline _dup(75) "_"
local title " Antimodes in circular KDE"
count if `deg' !=.
noi di as txt "`title', bw = " as res `h' as txt ", Ker = " as res `kc' as txt ", deg points = " as res r(N)
*noi di as txt "`title', bw = `hv', M = `mv', Ker = `kc'"
noi di as txt _dup(75) "-"
while `i'<`nuamo'+1 {
noi di as txt " Antimode ( " as res %4.0f `i' as txt " ) = " as res %12.4f `amodes'[`i']
local i = `i'+1
}
noi di as txt _dup(75) "_"
sort `deg'
}
}
qui if "`gen'" != "" {
restore, not
parse "`gen'", parse(" ")
gen `1' = `f' if `deg' != .
label var `1' ///
"Circular `klab' kernel density estimate, half-width `hlab'`=char(176)'"
gen `2' = `deg' if `deg' != .
label var `2' "Degrees `=char(176)'"
}
end