program define gbgfit_ll version 8.2 args lnf a b p q quietly replace `lnf' = ($S_mln)*ln(ibeta(`p',`q',(${S_mlz2}/`b')^`a'/(1+(${S_mlz2}/`b')^`a'))-ibeta(`p',`q',(${S_mlz1}/`b')^`a'/(1+(${S_mlz1}/`b')^`a'))) quietly replace `lnf' = ($S_mln)*ln(ibeta(`p',`q',(${S_mlz2}/`b')^`a'/(1+(${S_mlz2}/`b')^`a'))) if ${S_mlz1}==0 quietly replace `lnf' = ($S_mln)*ln(ibeta(`p',`q',(${S_mlz1}/`b')^`a'/(1+(${S_mlz1}/`b')^`a'))) if ${S_mlz2}==. tempvar sum qui g double `sum' =lngamma(sum(1+$S_mln)) qui count if `lnf'<. quietly replace `lnf' = `sum'[_N]/r(N) +`lnf' end exit