program define ci_marg_mu *! version 1.3 SRH 18 Oct 2008 version 10.0 syntax newvarlist(min=2 max=2) [if] [in] [, DOTS Level(int $S_level) Reps(int 1000)] if "`e(cmd)'" != "gllamm" { di in red "gllamm was not the last command" error 301 } local n200 = `reps'/200 if floor(`n200')<`reps'/200{ disp in re "reps must be a multiple of 200" exit 198 } tempname p b v bc vc T a cholv tempvar idno gen long `idno' = _n /* number of observation in tails */ local extreme = `n200'*(200-2*`level') timer clear 6 timer on 6 /* deal with if and in */ tempvar touse mark `touse' `if' `in' preserve qui keep if `touse' /* deal with estimates and constraints*/ matrix `b'=e(b) matrix `v'=e(V) if `e(const)'==1{ matrix `T'=e(T) matrix `a'=e(a) } else{ matrix `T'=I(colsof(`b')) matrix `a'=J(1,colsof(`b'),0) } matrix `vc' = `T''*`v'*`T' matrix `cholv' = cholesky(`vc') matrix `bc' = (`b'-`a')*`T' local vlist forvalues i=1/`extreme' { tempvar pm`i' local vlist `vlist' `pm`i'' mata: mf_sample_p() qui gllapred `pm`i'' if `touse', mu marg fsample from(`p') if "`dots'"!="" noi di in gr "." _c } if "`dots'"!="" noi disp " " timer off 6 quietly timer list 6 if "`dots'"!="" noi disp in gr r(t6) " seconds = " r(t6)/60 " minutes = " r(t6)/3600 " hours" if "`dots'"!="" noi disp " " timer clear 6 timer on 6 local extremep = `extreme' + 1 forvalues j=`extremep'/`reps'{ tempvar new mata: mf_sample_p() qui gllapred `new' if `touse', mu marg fsample from(`p') if "`dots'"!="" noi di in gr "." _c mata: rowshuffle(`extreme',"`touse'","`vlist' `new'") drop `new' } if "`dots'"!="" noi disp " " timer off 6 quietly timer list 6 if "`dots'"!="" noi disp in gr r(t6) " seconds = " r(t6)/60 " minutes = " r(t6)/3600 " hours" tokenize `varlist' local y "`1'" local lhalf = `extreme'/2 local uhalf = `lhalf' + 1 rename `pm`lhalf'' `1' rename `pm`uhalf'' `2' local lp = round((100-`level')/2,.1) local up = round(100-`lp',.1) label var `1' "`lp' percentile for gllapred, mu marg " label var `2' "`up' percentile for gllapred, mu marg " /* restore data */ tempfile file qui keep `idno' `1' `2' sort `idno' qui save "`file'", replace restore tempvar mrge sort `idno' qui merge `idno' using "`file'", _merge(`mrge') end version 10.0 mata: void rowshuffle(numeric scalar num, string scalar touse, string scalar varlist) { veclabs = tokens(varlist) // veclabs st_view(V,.,veclabs,touse) // sort rows and delete 26th element uptoi = num/2 fromi = uptoi + 2 nump = num + 1 for (i=1;i<=rows(V);i++){ newrow=sort(V[i,]',1)' newrow=newrow[1..uptoi],newrow[fromi..nump] V[i,1..num]=newrow } } void mf_sample_p() { cholV = st_matrix(st_local("cholv")) p = st_matrix(st_local("bc")) + (cholV * rnormal(cols(cholV), 1, 0, 1))' p = p*st_matrix(st_local("T"))' + st_matrix(st_local("a")) st_matrix(st_local("p"),p ) } end