*! version 1.0 93/07/22                 STB-18: sg23
program define bhatmesd
 version 3.0
 local varlist "req ex max(3)"
 local in "opt"
 if "`3'"=="" {
    di "syntax: bhatmesd freqvar diflogvar midpoivar [in range]" 
    exit}
 parse "`*'"
 quietly {
         tempvar touse
         gen byte `touse'=0
         replace `touse'=1 `in'
         regress `2' `3' if `touse'
         tempvar width mean sd b gauprob sumgaup sumsmfr en gaucom
         gen `width'=`3'[2]-`3'[1]
         gen `mean'=-1*_b[_cons]/_b[`3']+`width'/2
         gen `sd'=sqrt(`width'/(-1*_b[`3']))
         }
 gen `b'=_b[`3']
 if `b'>0 {
    di "The slope of selected points is >0" 
    exit}
 else 
 di "R-square = " %5.4f _result(7) _skip(10) "Adj R-square = " %5.4f _result(8) 
 di "Mean = " %5.4f `mean'
 di "s.d. = " %5.4f `sd'
 quietly {
 gen `gauprob'=`width'*(1/(`sd'*sqrt(2*_pi)))*exp(-.5*((`3'-`mean')/`sd')^2)
 egen `sumgaup'=sum(`gauprob') if `touse'
 egen `sumsmfr'=sum(`1') if `touse'
 gen `en'=`sumsmfr'[2]/`sumgaup'[2]
 replace `en'=round(`en',1)
 noisily di "component size = " `en'
 gen `gaucom'=`en'*`gauprob'
 sort `3'
 }
 graph `1' `gaucom' `3', xlab ylab c(ls) s(o.)
end