*! IHSU 1.1.0 04 July 2012. Based on the program by: * 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 12 May 2011 by Isaias Hazarmabeth Salgado-Ugarte * Revised 04 July 2012 by Isaias Hazarmabeth Salgado-Ugarte * This program calculates kernel density estimator of a series of * values according to the weight functions coded as follows: * 1 = Uniform * 2 = Triangle * 3 = Epanechnikov * 4 = Quartic (Biweight) * 5 = Triweight * 6 = Gaussian * 7 = Cosinus program define circkde version 5.0 local varlist "max(1)" local if "opt" local in "opt" local options "H(real 30) Kc(integer 6) NUMOdes MOdes NUAMOdes AMOdes noGraph * GENPDF(str) GENDEG(str)" parse "`*'" tempvar touse deg d e f qui { mark `touse' `if' `in' markout `touse' `varlist' count if `touse' local n = _result(1) gen `deg' = (_n / _N) * 360 gen `d' = . gen `e' = . gen `f' = . label var `deg' "Angle in degrees" 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 { replace `e'=cond(`e' >=`h',0,(_pi/4)*cos((_pi/2)*`e'/`h')) } *replace `e' = cond(`e' >= `h', 0, (1 - `e'^2 / `h'^2)^2) su `e', meanonly replace `f' = _result(18) / (`n' * `h') in `i' local i = `i' + 1 } if "`graph'" ~= "nograph" { if "`t1title'" ==""{ local t1title "Circular Kernel Density Estimation" local t1title "`t1title', bw = `h', k = `kc'" } graph `f' `deg', `options' t1("`t1title'") } if "`genpdf'" != "" { confirm new variable `genpdf' gen `genpdf' = `f' } if "`numodes'"~="" { tempvar difvar inmo sumo gen `difvar'=`f'[_n+1] - `f'[_n] gen `inmo' = 0 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" 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] gen `inamo' = 0 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 density estimation" 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' } } if "`gendeg'" != "" { confirm new variable `gendeg' gen `gendeg' = `deg' } end