*! version 1.3.0 14May2012 MLB program define marg_nb2cumul_prep, rclass version 10.1 syntax [if] , /// [ /// simvars(string) /// pi(string) /// simopts(string) /// noPArsamp /// ] /// ptheor(string) /// ytheor(string) /// pobs(varname) /// sims(integer) /// n(integer) /// y(varname) /// e(real) marksample touse tempvar lambda first fw alpha xg qui predict double `lambda' if `touse', xb qui replace `lambda' = exp(`lambda') qui predict double `alpha' if `touse', xb eq(#2) qui replace `alpha' = exp(`alpha') qui gen `xg' = . qui gen `ptheor' = . sum `y' if `touse', meanonly local min = r(min) local max = r(max) if "`parsamp'" == "" { tempname bname vname orig matrix `bname' = e(b) matrix `vname' = e(V) est store `orig' tempvar lambda2 alpha2 xg2 qui gen `lambda2' = . qui gen `alpha2' = . qui gen `xg2' = . forvalues i = 1/`sims' { local nvar : word `i' of `simvars' local pvar : word `i' of `pi' drop `lambda2' `alpha2' `xg2' marg_parsamp, bname(`bname') v(`vname') qui predict double `lambda2' if `touse', xb qui replace `lambda2' = exp(`lambda2') qui predict double `alpha2' if `touse', xb eq(#2) qui replace `alpha2' = exp(`alpha2') qui gen double `xg2' = rgamma(1/`alpha2', `alpha2') qui gen `nvar' = rpoisson(`lambda2'*`xg2') if `touse' sum `nvar' , meanonly local min = min(`min', `r(min)') local max = max(`max', `r(max)') sort `nvar' qui gen `pvar' = sum(`touse') qui replace `pvar' = cond(`touse',`pvar'/(`pvar'[_N]+1),.) quietly { by `nvar' : replace `pvar' = `pvar'[_N] } local gr "`gr' || scatter `pvar' `nvar', msymbol(oh) color(gs10) sort `simopts'" } qui est restore `orig' } else { forvalues i = 1/`sims' { drop `xg' local nvar : word `i' of `simvars' local pvar : word `i' of `pi' qui gen double `xg' = rgamma(1/`alpha', `alpha') qui gen `nvar' = rpoisson(`lambda'*`xg') if `touse' sum `nvar' , meanonly local min = min(`min', `r(min)') local max = max(`max', `r(max)') sort `nvar' qui gen `pvar' = sum(`touse') qui replace `pvar' = cond(`touse',`pvar'/(`pvar'[_N]+1),.) quietly { by `nvar': replace `pvar' = `pvar'[_N] } local gr "`gr' || scatter `pvar' `nvar', msymbol(oh) color(gs10) sort `simopts'" } } sort `lambda' `alpha' by `lambda' `alpha': gen byte `first' = cond(`touse', _n == 1, 0) qui by `lambda' `alpha' : gen float `fw' = _N if `touse' if _N < ( `max' + 1 ) { local n_old = _N qui set obs `= r(N) + 1 ' } qui gen `ytheor' = _n - 1 if _n <= ( `max' + 1 ) tempvar touse2 gen byte `touse2' = !missing(`ytheor') mata marg_nb2pp("`ytheor'", "`ptheor'", "`lambda'", "`alpha'", "`first'", "`fw'", "`touse2'") return local gr `"`gr'"' end