*!version 1.1 1999 Joseph Hilbe * first version 1.0.0 1993 Joseph Hilbe, Walter Linde-Zwirble (sg44: STB-28) * Beta-binomial random deviate generator * Example: rndbb 1000 200 0.2 0.05 * set obs 1000 * 200 is the denominator * 0.2 is p=a1/(a1+a2) * 0.05 is k, the dispersion factor k=1/(a1+a2+1) * This generator will return beta-binomial random deviates within the * following constraints. Although k can take any value from 0 to 1, in this * program we limit k because of the volotility of the distribution outside * this range of k. k must be as follows: k< the lessor of p' and (1-p')/2 * where p'=p if p<0.5, else p'=1-p. This should work well within the * physical representation of an overdispersed binomial. program define rndbb version 4.0 set type double cap drop xbb qui { local cases `1' set obs `cases' mac shift local m `1' mac shift local pp `1' mac shift local k `1' if `pp' < 0.5 { local p=`pp' } else { local p=1-`pp' } local eta=(1/`k')-1 if (`k' > `p') | (`k' > (1-`p')/2) { noi di in bl " k too large from program..." error } local varr=(`m'*`p'*(1-`p'))*(1+`k'*(`m'-1)) if `p'>0.4 { local peak=`m'*`p' local v=4*`varr' } else { local peak=`m'*`p'*(1-1.5*((`k'*`k')/(`p'*`p'))) local v=3*`varr' } #delimit ; local scale=exp(lngamma(`p'*`eta'+`peak')+ lngamma(`m'-`peak'+(1-`p') * `eta') - lngamma(`peak'+1) - lngamma(`m'-`peak'+1)); #delimit cr local sq=sqrt(`v') local am=`peak' tempvar ran1 ran2 ds ts sum1 t em y gen `ran1'=uniform() gen `ran2'=uniform() gen `ds'=1 gen `ts'=1 gen `em'=-1 gen `t'=-1 gen `y'=-1 egen `sum1'=sum(`ds') noi di in gr "( Generating " _c while `sum1' > 0 { replace `y' = sin(_pi*`ran1')/cos(_pi*`ran1') replace `em'= int(`sq'*`y'+`am') if `ds'==1 replace `ts'=0 if (( (0 > `em') | (`em'>`m') ) & (`ds'==1)) #delimit ; replace `t' = 0.66 * (1+`y'*`y')*exp(lngamma(`p'*`eta'+`em') + lngamma(`m'-`em'+(1-`p')*`eta')-lngamma(`em'+1) - lngamma(`m'-`em'+1))/`scale' if ((`ds'==1) & (`ts'==1)); #delimit cr replace `ds'=0 if ((`ran2'<`t') & (`ds'==1) & (`ts'==1)) replace `ran1'=uniform() replace `ran2'=uniform() replace `t'=-1 replace `ts'=1 drop `sum1' egen `sum1'=sum(`ds') noi di in gr "." _c } gen xbb = `em' replace xbb = `m'-`em' if (`pp'>`p') noi di in gr " )" noi di in bl "Variable " in ye "xbb " in bl "created." lab var xbb "Beta binomial random variable" set type float } end