I have been facing the same problem. There is an option that allows you to
specify the values at which you will compute the density.
The problem I was facing was to conduct an equivalence evaluation by
comparing the density functions I have deleted much material from the
program below in order to make the code a little simpler. This was for a 3
group problem in which I was testing that the groups weren't equal versus an
alternative of equality.
/* mn is the minimum value in the three groups, mx is the
maximum. */
scalar mn=mn1*(mn1<=mn2 & mn1<=mn3)+mn2*(mn2<=mn1 &
mn2<=mn3)+mn3*(mn3<=mn1 & mn3<=mn2)
scalar mx=mx1*(mx1>=mx2 & mx1>=mx3)+mx2*(mx2>=mx1 &
mx2>=mx3)+mx3*(mx3>=mx1 & mx3>=mx2)
/* Scale so that we integrate from less than the minimum to greater
than the maximum - this allows for negative or positive minimum or maximum
*/
scalar mn=1.25*mn*(mn<0)+.5*mn*(mn>0)
scalar mx=1.25*mx*(mx>0)+.5*mx*(mx<0)
/* mn is the minimum value of the observations - choose h to be the
increment for the estimation */
scalar h=(mx-mn)/(_N-1)
gen xx=mn+h*(_n-1)
/* Now find densities, normalize, and find area under the maximum */
kdensity y1,gen(xyz1 f1) at(xx) nograph
kdensity y2,gen(xyz2 f2) at(xx) nograph
kdensity y3,gen(xyz3 f3) at(xx) nograph
integ f1 xx,gen(Sf1) replace
/* this is to ensure that the densities integrate to 1 - it's a small point
*/
replace f1=f1/Sf1[_N]
integ f2 xx,gen(Sf2) replace
replace f2=f2/Sf2[_N]
integ f3 xx,gen(Sf3) replace
replace f3=f3/Sf3[_N]
/* finds maximum density */
egen fmax=rmax(f1 f2 f3)
integ fmax xx,gen(Sfmx) replace
/* this is compared to some critical value - I have used normal densities to
get some critical values. However, I think it will be better to find the
critical values using a bootstrap method. For a null of equality, it should
be much easier */
Peter A. Lachenbruch, Ph. D.
Director, Division of Biostatistics
OBE/CBER/FDA
Phone (301) 827-3320
FAX (301) 827-5218
*
* For searches and help try:
* http://www.stata.com/support/faqs/res/findit.html
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/
