*! 1.0.0 NJC & SPJ 13 Dec 2006
* Fit two-parameter inverse Gaussian distribution by ML

/*------------------------------------------------ playback request */
program invgaussfit, eclass byable(onecall)
	version 8.1
	if replay() {
		if "`e(cmd)'" != "invgaussfit" {
			di as err "results for invgaussfit not found"
			exit 301
		}
		if _by() error 190
		Display `0'
		exit `rc'
	}
	if _by() by `_byvars'`_byrc0': Estimate `0'
	else Estimate `0'
end

/*------------------------------------------------ estimation */
program Estimate, eclass byable(recall)
	syntax varlist(max=1) [if] [in] [fw aw] [,  ///
		MUvar(varlist numeric) LAMBDAvar(varlist numeric) ///
		Robust Cluster(varname) Level(integer $S_level) noLOG * ]

	local y "`varlist'"
	marksample touse
	markout `touse' `varlist' `muvar' `lambdavar' `cluster'

	qui count if `y' <= 0 & `touse'
	if r(N) {
		noi di " "
		noi di as txt "warning: {res:`y'} has `r(N)' values <= 0;" _c
		noi di as txt " not used in calculations"
	}
	qui replace `touse' = 0 if `y' <= 0

	qui count if `touse'
	if r(N) == 0 error 2000

	local title "ML fit of two-parameter inverse Gaussian dist'n"
	local nmu : word count `muvar'
	local nlambda : word count `lambdavar'

	local wtype `weight'
	local wtexp `"`exp'"'
	if "`weight'" != "" local wgt `"[`weight'`exp']"'

	if "`cluster'" != "" {
		local robust "robust"
		local clopt "cluster(`cluster')"
	}

	if "`level'" != "" local level "level(`level')"
        local log = cond("`log'" == "", "noisily", "quietly")

	mlopts mlopts, `options'
	global S_MLy "`y'"
	`log' ml model lf invgaussfit_lf (mu: `muvar') (lambda: `lambdavar') ///
		`wgt' if `touse' , maximize 				  ///
		collinear title(`title') `robust'       		  ///
		search(on) `clopt' `level' `mlopts' `stdopts' `modopts'

	eret local cmd "invgaussfit"
	eret local depvar "`y'"

	tempname b bmu blambda
	mat `b' = e(b)
	mat `bmu' = `b'[1,"mu:"]
	mat `blambda' = `b'[1,"lambda:"]
	eret matrix b_mu = `bmu'
	eret matrix b_lambda = `blambda'
	eret scalar length_b_mu = 1 + `nmu'
	eret scalar length_b_lambda = 1 + `nlambda'

	if "`muvar'`lambdavar'" == ""  {
		tempname e
		mat `e' = e(b)
		eret scalar mu = `e'[1,1]
		eret scalar lambda = `e'[1,2]
	}

	Display, `level' `diopts'
end

program Display
	syntax [, Level(int $S_level) *]
	local diopts "`options'"
	ml display, level(`level') `diopts'
	if `level' < 10 | `level' > 99 local level = 95
end

*! NJC 1.0.0 7 Dec 2006
program invgaussfit_lf
	version 8.1
	args lnf mu lambda
	qui replace `lnf' = ///
	0.5 * ln(`lambda') - ln(2 * _pi) - 3 * ln($S_MLy) ///
	- (`lambda'/(2 * `mu'^2 * $S_MLy) * ($S_MLy - `mu')^2)
end