*! version 2.1.0 12feb2012 *! author: Partha Deb * version 1.0.1 24apr2007 * version 1.0.0 06mar2007 program fmm_poisson_lf version 9.2 forvalues i = 1/\$fmm_components { local L_xb `"`L_xb' xb`i'"' local L_exb `"`L_exb' exb`i'"' local L_fxb `"`L_fxb' fxb`i'"' local L_pr `"`L_pr' pr`i'"' local L_gb `"`L_gb' gb`i'"' } forvalues i = 1/`=2*\$fmm_components-1' { local L_gL `"`L_gL' gL`i'"' forvalues j = `i'/`=2*\$fmm_components-1' { local L_h `"`L_h' h`i'`j'"' local L_nh `"`L_nh' nh`i'`j'"' } } forvalues i = 1/`=\$fmm_components-1' { local L_lpr `"`L_lpr' lpr`i'"' forvalues j = 1/`=\$fmm_components-1' { local L_ga `"`L_ga' ga`i'`j'"' } } // model arguments and temporary variables args todo b lnf g negH `L_gL' tempname `L_xb' `L_exb' `L_fxb' `L_lpr' `L_pr' /// den prob gi `L_gb' `L_ga' hij `L_h' `L_nh' // set up equations forvalues i=1/\$fmm_components { mleval `xb`i'' = `b', eq(`i') qui gen double `exb`i'' = exp(`xb`i'') } qui gen double `den' = 1 forvalues i=1/`=\$fmm_components-1' { mleval `lpr`i'' = `b', eq(`=\$fmm_components+`i'') qui replace `den' = `den' + exp(`lpr`i'') } forvalues i=1/`=\$fmm_components-1' { qui gen double `pr`i'' = exp(`lpr`i'')/`den' } qui gen double `pr\$fmm_components' = 1 forvalues i=1/`=\$fmm_components-1' { qui replace `pr\$fmm_components' = `pr\$fmm_components' - `pr`i'' } // calculate the likelihood function qui gen double `prob' = 0 forvalues i=1/\$fmm_components { qui gen double `fxb`i'' = exp(\$ML_y*`xb`i'' - `exb`i'' /// - lngamma(\$ML_y+1)) qui replace `prob' = `prob' + `pr`i''*`fxb`i'' } mlsum `lnf' = ln(`prob') // CALCULATE GRADIENT TERMS // gradient bi forvalues i = 1/\$fmm_components { qui gen double `gb`i'' = \$ML_y-`exb`i'' // density specific qui replace `gL`i'' = (`pr`i'' * `fxb`i'' * `gb`i'')/`prob' } // gradient prj forvalues j = 1/`=\$fmm_components-1' { local m = `=\$fmm_components+`j'' qui replace `gL`m'' = 0 forvalues k = 1/`=\$fmm_components-1' { qui gen double `ga`k'`j'' = ((`j'==`k') - `pr`j'') qui replace `gL`m'' = `gL`m'' + `pr`k''*`ga`k'`j'' /// *(`fxb`k'' - `fxb\$fmm_components')/`prob' } } // collect gradient terms into vector local np = colsof(`b') local c 1 matrix `g' = J(1,`np',0) forvalues i = 1/`=2*\$fmm_components-1' { mlvecsum `lnf' `gi' = `gL`i'', eq(`i') matrix `g'[1,`c'] = `gi' local c = `c' + colsof(`gi') } // CALCULATE HESSIAN TERMS // hessian - b terms local c 1 qui gen double `hij' = . forvalues i = 1/\$fmm_components { // hessian (bi,bi) qui replace `hij' = -`exb`i'' // density specific qui gen double `h`i'`i'' = `pr`i''/`prob'*`fxb`i'' /// *(-`gL`i''*`gb`i'' + `gb`i''^2 + `hij') mlmatsum `lnf' `nh`i'`i'' = -`h`i'`i'', eq(`i',`i') // hessian (bi,bj) if (`i'<\$fmm_components) { forvalues j = `=`i'+1'/\$fmm_components { qui gen double `h`i'`j'' = `pr`i''/`prob'*(-`gL`j''*`fxb`i''*`gb`i'') mlmatsum `lnf' `nh`i'`j'' = -`h`i'`j'', eq(`i',`j') } } // hessian (bi,prj) if (`i'<\$fmm_components) { forvalues j = 1/`=\$fmm_components-1' { local m = `=\$fmm_components+`j'' qui gen double `h`i'`m'' = -`gL`m''*`gL`i'' /// + 1/`prob'*`fxb`i''*`gb`i''*`pr`i''*`ga`i'`j'' mlmatsum `lnf' `nh`i'`m'' = -`h`i'`m'', eq(`i',`m') } } else { forvalues j = 1/`=\$fmm_components-1' { local m = `=\$fmm_components+`j'' qui gen double `h`i'`m'' = -`gL`m''*`gL`i'' forvalues k = 1/`=\$fmm_components-1' { qui replace `h`i'`m'' = `h`i'`m'' /// - 1/`prob'*`fxb`i''*`gb`i''*`pr`k''*`ga`k'`j'' } mlmatsum `lnf' `nh`i'`m'' = -`h`i'`m'', eq(`i',`m') } } } // hessian - pr terms // hessian (prj,pri) forvalues i = 1/`=\$fmm_components-1' { forvalues j = `i'/`=\$fmm_components-1' { local m = `=\$fmm_components+`i'' local n = `=\$fmm_components+`j'' qui gen double `h`m'`n'' = -`gL`m''*`gL`n'' qui replace `hij' = -`pr`i''*((`i'==`j') - `pr`j'') forvalues k = 1/`=\$fmm_components-1' { qui replace `h`m'`n'' = `h`m'`n'' /// + 1/`prob'*`pr`k''*(`fxb`k'' - `fxb\$fmm_components') /// *(`ga`k'`i''*`ga`k'`j'' + `hij') } mlmatsum `lnf' `nh`m'`n'' = -`h`m'`n'', eq(`m',`n') } } // collect hessian terms into matrix local np = colsof(`b') local r 1 matrix `negH' = J(`np',`np',0) forvalues i = 1/`=2*\$fmm_components-1' { local c = `r' forvalues j = `i'/`=2*\$fmm_components-1' { matrix `negH'[`r',`c'] = `nh`i'`j'' if (`j'>`i') { matrix `negH'[`c',`r'] = `nh`i'`j''' } local c = `c' + colsof(`nh`i'`j'') } local r = `r' + rowsof(`nh`i'`i'') } end