*! version 2.1.0 12feb2012
*! author: Partha Deb
* version 1.0.1  24apr2007
* version 1.0.0 06mar2007

program fmm_gamma_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_lalp `"`L_lalp' lnalpha`i'"'
		local L_alp `"`L_alp' alpha`i'"'
		local L_gb `"`L_gb' gb`i'"'
		local L_gs `"`L_gs' gs`i'"'
	}

	forvalues i = 1/`=3*$fmm_components-1' {
		local L_gL `"`L_gL' gL`i'"'
		forvalues j = `i'/`=3*$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_lalp' `L_alp' `L_lpr' `L_pr' ///
						den prob gi `L_gb' `L_ga' `L_gs' 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'')
		mleval `lnalpha`i'' = `b', eq(`=2*$fmm_components-1+`i'')
		qui gen double `alpha`i'' = exp(`lnalpha`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'' = gammaden(`alpha`i'',`exb`i'',0,$ML_y)
		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'' = -`alpha`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'
		}
	}

	// gradient alphai
	forvalues i=1/$fmm_components {
		local k = `=2*$fmm_components-1+`i''
		qui gen double `gs`i'' = (-digamma(`alpha`i'') - `xb`i'' ///
															+ log($ML_y))*`alpha`i''	// density specific
		qui replace `gL`k'' = (`pr`i'' * `fxb`i'' * `gs`i'')/`prob'
	}

	// collect gradient terms into vector
	local np = colsof(`b')
	local c 1
	matrix `g' = J(1,`np',0)

	forvalues i = 1/`=3*$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' = -$ML_y/`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 (bi,alpha)
		forvalues j = 1/$fmm_components {
			local k = `=2*$fmm_components-1+`j''
			// hessian (bi,alphai)
			if (`i'==`j') {
				qui replace `hij' = -`alpha`j''			// density specific
				qui gen double `h`i'`k'' = `pr`i''/`prob'*`fxb`i'' ///
					*(-`gL`k''*`gb`i'' + `gs`j''*`gb`i'' + `hij')
				mlmatsum `lnf' `nh`i'`k'' = -`h`i'`k'', eq(`i',`k')
			}
			else {
			// hessian (bi,alphaj)
				qui gen double `h`i'`k'' = `pr`i''/`prob'*(-`gL`k''*`fxb`i''*`gb`i'')
				mlmatsum `lnf' `nh`i'`k'' = -`h`i'`k'', eq(`i',`k')
			}
		}
	}

	// 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')
		}
	}

	// hessian - alpha terms
	forvalues i = 1/$fmm_components {
		// hessian w.r.t. alphai
		local m = `=2*$fmm_components-1+`i''
		qui replace `hij' = (-digamma(`alpha`i'') - `xb`i'' + log($ML_y) ///
								- trigamma(`alpha`i'')*`alpha`i'')*`alpha`i''	// density specific
		qui gen double `h`m'`m'' = `pr`i''/`prob'*`fxb`i'' ///
			*(-`gL`m''*`gs`i'' + `gs`i''^2 + `hij')
		mlmatsum `lnf' `nh`m'`m'' = -`h`m'`m'', eq(`m',`m')
		// hessian w.r.t. alphaj (cross partials)
		forvalues j = `=`i'+1'/$fmm_components {
			local n = `=2*$fmm_components-1+`j''
			qui gen double `h`m'`n'' = `pr`i''/`prob'*(-`gL`n''*`fxb`i''*`gs`i'')
			mlmatsum `lnf' `nh`m'`n'' = -`h`m'`n'', eq(`m',`n')
		}

	// hessian (alphai,prj)
		if (`i'<$fmm_components) {
			forvalues j = 1/`=$fmm_components-1' {
				local n = `=2*$fmm_components-1+`i''
				local m = `=$fmm_components+`j''
				qui gen double `h`m'`n'' = -`gL`m''*`gL`n'' ///
						+ 1/`prob'*`fxb`i''*`gs`i''*`pr`i''*`ga`i'`j''
				mlmatsum `lnf' `nh`m'`n'' = -`h`m'`n'', eq(`m',`n')
			}
		}
		else {
			forvalues j = 1/`=$fmm_components-1' {
				local n = `=2*$fmm_components-1+`i''
				local m = `=$fmm_components+`j''
				qui gen double `h`m'`n'' = -`gL`m''*`gL`n''
					forvalues k = 1/`=$fmm_components-1' {
						qui replace `h`m'`n'' = `h`m'`n'' ///
							- 1/`prob'*`fxb`i''*`gs`i''*`pr`k''*`ga`k'`j''
					}
				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/`=3*$fmm_components-1' {
		local c = `r'
		forvalues j = `i'/`=3*$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