*! version 3.4.1 18Mar97. STB-36 sbe13.2
program define xriml
version 4.0
if "`*'"=="" | substr("`1'",1,1)=="," { 
	mx_out mlreg	/* Note: local version of mx_out (ml mlout) */
	exit
}
loc varlist "req ex min(1) max(2)"
loc if "opt"
loc in "opt"
#delimit ;
local hidopts "EXpx(str) ORIgin(str) Zero noTIdy TRace ZCI PARam(str)";
loc options "FP(str) COVArs(str) CV noGRaph noLEAve 
 CEntile(str) SE DIst(str) noSCALing LTolerance(real .001) 
 INit(str) SAVing(str) `hidopts'";
#delimit cr
loc small 1e-6
parse "`*'"
if "`scaling'"!="noscaling" { local showtsf "noi" }
else local showtsf "qui"
if "`dist'"=="" {
	di in red "invalid distribution"
	exit 198
}
if "`centile'"=="" { loc centile "3 97" }
parse "`centile'", parse(" ,")
local centile
loc i 0
while "`1'"!="" {
	if "`1'"!="," {
		loc nc=`nc'+1
		loc c`nc' `1'
		loc C "`C' __C`1'"
		loc centile "`centile' `1'"
	}
	mac shift
}
loc dist = lower(substr("`dist'",1,2))
if "`dist'"=="n"|"`dist'"=="no" {
	loc dnum 0
	loc dist n
}
else {
	loc d= 1*("`dist'"=="sl")+2*("`dist'"=="pn")+3*("`dist'"=="en") /*
	 */   +7*("`dist'"=="mp")+8*("`dist'"=="me")
	if `d'==0 {
		di in red "invalid distribution `dist'"
		exit 198
	}
	else loc dnum `d'
}
loc pars3=(`dnum'>=3)
loc pars4=(`dnum'>=4)
if `dnum'>0 {
	loc gamma "g"
}
if `pars4' {
	loc delta "d"
}
loc dupper=upper("`dist'")
loc dt0 "Normal"
loc dt1 "Shifted Log-Normal"
loc dt2 "Power-Normal"
loc dt3 "Exponential-Normal"
loc dt7 "Modulus power-Normal"
loc dt8 "Modulus exp-Normal"
if "`expx'"!="" { loc expx "expx(`expx')" }
if "`origin'"!="" { loc origin "origin(`origin')" }
if "`detail'"!="" { loc detail "noisily" }
if "`cv'"!="" { loc cv "cv" }

if "`param'"!="" {
	parse "`param'", parse(" ")
	while "`1'"!="" {
		loc fl=lower(substr("`1'",1,2))
		if "`fl'"=="cv" 	{ loc cv "cv" }
		else if "`fl'"=="ln" 	{ loc lns "lns" }
		else if "`fl'"=="ta" 	{ loc tau "tau" }
		else if "`fl'"=="tr" 	{ loc trunc "trunc" } /*!! not SU */
		else {
			di in red "invalid param() option `1'"
			exit 198
		}
		mac shift
	}
}
/*
	Extract base variables for each parameter
*/
_jprxrpa "`covars'" ms`gamma'`delta' "covariates"
loc ncl $S_1
loc i 1
while `i'<=`ncl' {
	loc curve ${S_p`i'}
	loc `curve'base "${S_c`i'}" /* basevars */
	loc i = `i'+1
}
/*
	Parameter initialisation: gamma (& delta) only at present
*/
_jprxrpa "`init'" `gamma'`delta' initializations
loc ncl $S_1 /* ncl is at most 2, but following code is general */
loc i 1
while `i'<=`ncl' {
	loc curve ${S_p`i'}
	loc init "${S_c`i'}"
/*
	!! ? Can be initialised with a var or a num
*/
	cap confirm var `init'
	if _rc { confirm num `init' }
	loc `curve'init `init'
	loc i = `i'+1
}
/*
	Default initial shape params depend on distribution
*/
if "`ginit'"=="" {
	loc ginit 0
	if `dnum'==0		{ /* normal */
		loc ginit 1
		loc gfix "gfix"
	}
	else if `dnum'==2	{ loc ginit 1 }    /* PN  */
	else if `dnum'==3	{ loc ginit 0.01 } /* EN  */

	else if `dnum'==7	{ loc ginit 1 }    /* MPN */
	else if `dnum'==8	{ loc ginit -0.2 } /* MEN */
}
if "`dinit'"=="" & `pars4' {
	if `dnum'==7 | `dnum'==8 { /* MPN, MEN */
		loc dinit 1
	}
	else loc dinit 0
}

parse "`varlist'", parse(" ")
tempvar y touse M Z
loc lhs `1'
loc rhs `2'
mark `touse' `if' `in'
markout `touse' `lhs' `rhs' `mbase' `sbase' `gbase' `dbase'
/*
	Check for null entries for M and S curves in `fp';
	if so, substitute default model specifications.

!!	Code commented out, default now constants for everything

_jprxrpa "`fp'" ms`gamma'`delta' "FP model items"
loc ncl $S_1 /* could be 0 */
loc i 1
loc mdef "m:powers 1 2,"
loc sdef "s:powers 1,"
while `i'<=`ncl' {
	loc curve ${S_p`i'}
	loc `curve'def
	loc i=`i'+1
}
loc fp "`mdef'`sdef'`fp'" /* does not include `gdef' and `ddef', but could */
*/
/*
	Extract FP model specifications.
*/
_jprxrpa "`fp'" ms`gamma'`delta' "FP model items"
loc ncl $S_1 /* could be 0 */
loc i 1
while `i'<=`ncl' {
	loc curve ${S_p`i'}
	loc spec "${S_c`i'}"
	parse "`spec'", parse(" ")
	loc first = lower("`1'")
	if "`first'"=="fix" {
		if "`curve'"=="m" | "`curve'"=="s" {
			di in red "invalid `spec'"
			exit 198
		}
		loc `curve'fix "`curve'fix"
		mac shift
		if "`*'"!="" {
			cap confirm var `*'
			if _rc { confirm num `*' }
			loc `curve'init `*'
		}
	}
	else {
		if "`rhs'"=="" {
			di in red "no xvar, cannot specify FP model"
			exit 198
		}
		if "`first'"=="powers" { mac shift }
		`showtsf' fracgen `rhs' `*', `scaling' /*
		 */ name(X`curve') replace `expx' `origin' `zero'
		loc `curve'vars $S_1
		loc `curve'fixpow `*'
	}
	loc i = `i'+1
}
/*
	Next 6 stmts set globals for use by _`dist'll.ado.
	Possibly could be made into _dta chars.
*/
global S_trunc = "`trunc'"!=""
global S_cv = "`cv'"!=""
global S_tau = "`tau'"!=""
global S_lns = "`lns'"!=""
global S_gamma `ginit'
global S_delta `dinit'
global S_dist `dist' /* for default use by centcalc and centcal2 */
global S_gfix = "`gfix'"!=""
global S_dfix = "`dfix'"!=""
quietly {
/*
	Calculate and set up initial parameter values for ML fitting
	!! needs thought about how/if to incorporate initialisation for M and S
*/
	tempname s g d mu b f V
/*
	Equation management. In Stata's eq.ado, the current highest equation
	number is stored in $S_eqnum. If you "eq drop" the most recently def-
	ined equation, the value of $S_eqnum is decremented by one, as hoped.
	If you drop something below the most recent, $S_eqnum is not reduced
	and problems eventually arise when $S_eqnum exceeds 99.

	So here we drop in the order Dcurve, Gcurve, Scurve, Mcurve,
	as Dcurve (if defined) would be the most recent from a previous run.
*/
	cap eq drop Dcurve
	cap eq drop Gcurve
	cap eq drop Scurve
	cap eq drop Mcurve
	eq Mcurve : `lhs' `mvars' `mbase'
	eq Scurve : `svars' `sbase'

	reg `lhs' `mvars' `mbase' if `touse', nohead
	matrix `mu' = get(_b)
	matrix colnames `mu' = Mcurve:
	loc obs = _result(1)
	if "`absres'"!="" {
		loc scale = _result(9) /* initialisation of sigma (_cons) */
		if "`cv'"!="" { /* CV parametrization */
			sum `lhs' if `touse'
			loc scale = `scale'/_result(3)
		}
		if "`lns'"!="" { loc scale = ln(`scale') }
		matrix `s' = (`scale')
		matrix colnames `s' = Scurve:_cons
	}
	else {
/*
	Initialize S via absolute residuals
*/
		predict `M' if `touse'
		if "`lns'"!="" {
/*
	-1.96... is digamma(0.5), Abramowitz & Stegun p258
*/
			scalar `f' = -0.5*(-1.963510026+ln(2))
			if "`cv'"!="" {
				gen `Z' = ln(abs(`lhs'/`M'-1))+`f'
			}
			else {
				gen `Z' = ln(abs(`lhs'-`M'))+`f'
			}
		}
		else {
			scalar `f' = sqrt(_pi/2)
			if "`cv'"!="" {
				gen `Z' = `f'*abs(`lhs'/`M'-1)
			}
			else {
				gen `Z' = `f'*abs(`lhs'-`M')
			}
		}
		reg `Z' `svars' `sbase', depname("Abs. res") nohead
		matrix `s' = get(_b)
		matrix colnames `s' = Scurve:
		drop `M' `Z'
	}
	matrix `mu' = `mu',`s'
	loc depvs "10"
	loc mods "Mcurve Scurve"
	if "`gfix'"=="" { /* gamma to be estimated */
		eq Gcurve : `gvars' `gbase'
		matrix `g' = (`ginit')
		matrix colnames `g' = Gcurve:_cons
		matrix `mu' = `mu',`g'
		loc depvs "`depvs'0"
		loc mods "`mods' Gcurve"
	}
	if `pars4' & ("`dfix'"=="") { /* delta to be estimated */
		eq Dcurve : `dvars' `dbase'
		matrix `d' = (`dinit')
		matrix colnames `d' = Dcurve:_cons
		matrix `mu' = `mu',`d'
		loc depvs "`depvs'0"
		loc mods "`mods' Dcurve"
	}
	ml begin
	ml function xri_`dist'll
	ml method lf
	ml model `b' = `mods', depv(`depvs') from(`mu')
	ml sample `touse' `if' `in'
	noi ml max `f' `V', ltol(`ltolera') `trace'
	tempname v
	matrix `v' = `V'
	ml post mlreg, title(`dt`dnum'' Regression)
	loc dev = -2*$S_E_ll	/* +`obs'*ln(2*_pi) now unnecessary */

	tempvar S G z
	predict `M' if `touse', equation(Mcurve)
	predict `S' if `touse', equation(Scurve)
	if "`gfix'"!="" {
		gen `G' = `ginit' if `touse'
	}
	else {
		predict `G' if `touse', equation(Gcurve)
	}
	if `pars4' {
		tempvar D
		if "`dfix'"!="" {
			gen `D' = `dinit' if `touse'
		}
		else {
			predict `D' if `touse', equation(Dcurve)
		}
	}
	if "`se'"!="" {
		tempvar se1 se2 se3
		predict `se1' if `touse', equation(Mcurve) stdp
		predict `se2' if `touse', equation(Scurve) stdp
		if "`gfix'"=="" {
			predict `se3' if `touse', equation(Gcurve) stdp
		}
		if `pars4' & ("`dfix'"=="") {
			tempvar se4
			predict `se4' if `touse', equation(Dcurve) stdp
		}
	}
/*
	Centiles (and standard errors)
*/
	if `pars4' {
		loc dd "delta(`D')"
	}
	centcalc `M' `S', centile(`centile') prefix(__C) /*
	 */ dist(`dist') gamma(`G') `dd' `cv' `lns' `tau' `trunc'
/*
	SEs of centile values.
*/
	loc no_x=("`rhs'`covars'"=="")
	if "`se'"!="" {
		loc i 0
		while `i'<`nc' {
		    loc i=`i'+1
		    tempvar c`i'se
		    gen `c`i'se'=.
		    noi di in gr _n "Calculating SE of `c`i''th centile ..."
		    if `pars4' {
			_se4 `dnum' `touse' /*
			 */ "`mvars'" "`svars'" "`gvars'" "`dvars'" /*
			 */ "`mbase'" "`sbase'" "`gbase'" "`dbase'" /*
			 */ `M' `S' `G' `D' `v' "`cv'" "`tau'" "`lns'" /*
			 */ "`gfix'" "`dfix'" `c`i'' `c`i'se' `no_x'
		    }
		    else {
			_se3 `dnum' `touse' "`mvars'" "`svars'" "`gvars'" /*
			 */ "`mbase'" "`sbase'" "`gbase'" /*
			 */ `M' `S' `G' `v' "`cv'" "`tau'" "`lns'" "`gfix'" /*
			 */ `c`i'' `c`i'se' `no_x'
		    }
		}
	}
/*
	Standardized residuals (Z-scores), calc by *ll.ado
*/
	gen `Z' = __U	/* convert __U to float */
	if "`zci'"!="" {
		tempvar yse zl zu
		gen `yse'=.
		if `pars4' {
			_se4 `dnum' `touse' /*
			 */ "`mvars'" "`svars'" "`gvars'" "`dvars'" /*
			 */ "`mbase'" "`sbase'" "`gbase'" "`dbase'" /*
			 */ `M' `S' `G' `D' `v' "`cv'" "`tau'" "`lns'" /*
			 */ "`gfix'" "`dfix'" `Z' `yse' `no_x'
		}
		else {
			_se3 `dnum' `touse' "`mvars'" "`svars'" "`gvars'" /*
			 */ "`mbase'" "`sbase'" "`gbase'" /*
			 */ `M' `S' `G' `v' "`cv'" "`tau'" "`lns'" "`gfix'" /*
			 */ `Z' `yse' `no_x'
		}
		loc z=-invnorm((100-$S_level)/200)
		gen `zl'=`lhs'-`z'*`yse'
		gen `zu'=`lhs'+`z'*`yse'
		global S_mldepn `zl'
		xri_`dist'll `yse' `M' `S' `G' `D' /* `yse' is a dummy here */
		replace `zl'=__U
		global S_mldepn `zu'
		xri_`dist'll `yse' `M' `S' `G' `D'
		replace `zu'=__U
	}
	drop __U
}
/*
	Graph
*/
if "`graph'"!="nograph" & "`rhs'"!="" {
	if "`saving'"!="" { local saving "saving(`saving')" }
	if "`mbase'`sbase'`gbase'`dbase'"!="" { loc cs "s(.pppppppp)" }
	else loc cs "c(.ssssssss) s(oiiiiiiii)"
	loc t1title "`dt`dnum'' model for `lhs'"
	loc t2title "`centile' centiles"
	graph `lhs' `M' `C' `rhs' if `touse', t1title("`t1title'") xlab ylab /*
	 */ t2title("`t2title'") l1("`lhs'") b2("`rhs'") sort `saving' `cs'
}
/*
	Save variables if required, and name/label them.
*/
if "`leave'"=="" {
	loc prog "_ml"
	loc vn Z`prog'
	cap drop `vn'
	rename `Z' `vn'
	lab var `vn' "`dupper' Z-scores"

	if "`cv'"!="" { loc cvl " [cv]" }
	if "`lns'"!="" { loc lnsl "ln " }

	loc a1 m
	loc a2 s

	loc b1 "location"
	loc b2 "`lnsl'scale`cvl'"

	if `pars3' {
		loc b3 "shape"
		loc a3 g
	}
	if `pars4' {
		loc b3 "shape [gamma]"
		loc a4 d
		loc b4 "shape [delta]"
	}
	loc i 1
	loc p=2+`pars3'+`pars4'
	while `i'<=`p' {
	    loc a `a`i''
	    loc A = upper("`a'")
	    loc b `b`i''
	    loc vn `A'`prog'
	    cap drop `vn'
	    rename ``A'' `vn'
	    if "``a'fix'"!=""  		 { loc add "fixed at ``a'init'" }
	    else if "``a'fixpow'"!="" 	 { loc add "powers ``a'fixpow'" }
	    else 			 { loc add "constant" }
	    lab var `vn' "`dupper' `b': `add'"
	    if "`se'"!="" & (`i'<=2 | (`i'==3 & "`gfix'"=="") /*
	     */ | ((`i'==4) & `pars4' & "`dfix'"=="")) {
		loc vn se`vn'
		cap drop `vn'
	    	rename `se`i'' `vn'
	    	lab var `vn' "`dupper' se[`b']"
	    }
	    loc i = `i'+1
	}
	loc i 0
	while `i'<`nc' {
		loc i=`i'+1
		loc cent `c`i''
		loc vn=substr("C`cent'`prog'",1,8)
		cap drop `vn'
		rename __C`cent' `vn'
		lab var `vn' "`dupper' `cent' centile"
		if "`se'"!="" {
			loc vn=substr("se`vn'",1,8)
			cap drop `vn'
			rename `c`i'se' `vn'
			lab var `vn' "`dupper' se[`cent' centile]"
		}
	}
	if "`zci'"!="" {
		loc vn "Zl`prog'"
		cap drop `vn'
		rename `zl' `vn'
		lab var `vn' "`dupper' Z lower ${S_level}% conf. limit"
		loc vn "Zu`prog'"
		cap drop `vn'
		rename `zu' `vn'
		lab var `vn' "`dupper' Z upper ${S_level}% conf. limit"
}
mx_out mlreg
di in gr _n "  Final deviance = " %9.3f in ye `dev' /*
 */ in gr " (" in ye `obs' in gr " observations.)"
if "`tidy'"!="notidy" {
	if "`mvars'"!="" { drop `mvars' }
	if "`svars'"!="" { drop `svars' }
	if "`gvars'"!="" { drop `gvars' }
	if "`dvars'"!="" { drop `dvars' }
}
global S_1 `dev'
end

program define _se3 /* SE for estimated centiles, 2 and 3 parameter models */
	* Updated cv case 17/3/97
	loc dnum  `1'	/* distribution number */
	loc touse `2'	/* filter for original obs */
	loc vars1 `3'	/* covariates for M-curve */
	loc vars2 `4'	/* covariates for S-curve */
	loc vars3 `5'	/* covariates for G-curve */
	loc base1 `6'	/* base vars  for M-curve */
	loc base2 `7'	/* base vars  for S-curve */
	loc base3 `8'	/* base vars  for G-curve */
	loc M `9'
	loc S `10'
	loc G `11'
	loc V `12'	/* VCE matrix of parameters */
	loc cv = "`13'"!=""
	loc tau = "`14'"!=""
	loc lns = "`15'"!=""
	loc gfix = "`16'"!=""
	loc centile `17'
	loc cse `18'	/* to hold SE, assumed to exist */
	loc no_x `19'
* Treat normal as boxcox

	if `dnum'==0 { loc dnum 2 }

* Put names of covariates into xn*

	loc small 1e-6

	loc p 1	/* parameter no.: 1=m, 2=s, 3=g. */
	loc npar 0	/* npar counts the covariates, including _cons's */
	while `p'<=3 {
		loc no_gpar = `gfix' & (`p'==3)
		if `no_gpar' { loc n3 0 }
		else loc n`p' 1
/*
	Put names of the FP and base covariates into xn*
*/
		if "`vars`p''`base`p''"!="" {
			parse "`vars`p'' `base`p''", parse(" ")
			while "`1'"!="" {
				loc n`p'=`n`p''+1
				loc npar=`npar'+1
				loc xn`npar' "`1'"
				mac shift
			}
		}
		if !`no_gpar' { loc npar=`npar'+1 }	/* add one for _cons */
		loc p=`p'+1
	}
/*
	Calc derivatives wrt params.
*/
	tempvar h1 h2 h3 A B
	cap confirm var `centile'
	if _rc {
		tempname uq
		scalar `uq'=invnorm(`centile'/100)
	}
	else local uq `centile' /* uq contains Z-scores (zci option) */
	if `lns' {	/* log S parameterization */
		tempvar s
		gen `s' = exp(`S')
	}
	else loc s `S'
	if `dnum'==1 { /* SLN */
	    gen `A' = exp(`G'*`uq')
	    gen `h1' = 1
	    if `cv' {
		gen `h2' = cond(abs(`G')<`small', `M'*`uq', /* 
		*/ `M'*(`A'-1)/`G' )
	    }
	    else gen `h2' = cond(abs(`G')<`small', `uq', (`A'-1)/`G' )
	    if !`gfix' {
		if `cv' {
			gen `h3' = cond(abs(`G')<`small', /*
			 */ 0.5*`s'*`M'*`uq'^2, /*
			 */ (`A'*(`uq'-1/`G')+1/`G')*`s'*`M'/`G' )
		}
		else {
			gen `h3' = cond(abs(`G')<`small', /*
			 */ 0.5*`s'*`uq'^2, /*
			 */ (`A'*(`uq'-1/`G')+1/`G')*`s'/`G' )
		}
	    }
	}
	else if `dnum'==2 {	/* PN */
	    if `tau' {		/* tau parameterization */
		tempvar t	/* lambda */
		if `cv' {
		    gen `t' = 1-`G'/`s'
		}
		else gen `t' = 1-`G'*`M'/`s'
	    }
	    else loc t `G'
	    if `cv' {
	    	gen `B' = `uq'*`s'
	    }
	    else gen `B' = `uq'*`s'/`M'
	    gen `A' = (1+`t'*`B')^(1/`t')
	    gen `h1' = cond(abs(`t')<`small', (1-`B')*exp(`B'), /*
		*/ `A'*(1-`B'/`A'^`t'))
	    if `cv' {
		gen `h2' = cond(abs(`t')<`small', `M'*`uq'*exp(`B'), /*
		*/ `M'*`uq'*`A'^(1-`t'))
	    }
	    else {
		gen `h2' = cond(abs(`t')<`small', `uq'*exp(`B'), /*
		*/ `uq'*`A'^(1-`t'))
	    }
	    if !`gfix' {
		gen `h3' = cond(abs(`t')<`small', -0.5*`M'*`B'^2*exp(`B'), /*
		 */ `M'*`A'*(`B'/`A'^`t'-log(`A'))/`t')
	    }
	}
	else if `dnum'==3 { /* EN */
	    gen `A' = 1+`G'*`uq'
	    gen `h1' = 1
	    if `cv' {
		    gen `h2' = cond(abs(`G')<`small', `M'*`uq', /*
		    */ `M'*log(`A')/`G' )
	    }
	    else gen `h2' = cond(abs(`G')<`small', `uq', log(`A')/`G' )
	    if !`gfix' {
		if `cv' {
			gen `h3' = cond(abs(`G')<`small', /*
			 */ -0.5*`s'*`M'*`uq'^2, /*
			 */ (1-1/`A'-log(`A'))*`s'*`M'/`G'^2 )
		}
		else {
			gen `h3' = cond(abs(`G')<`small', /*
			 */ -0.5*`s'*`uq'^2, /*
			 */ (1-1/`A'-log(`A'))*`s'/`G'^2 )
		}
	    }
	}
* Calc SE of centile, element-by-element
	tempname xx h hbeta hVh se
	matrix `hbeta' = J(1,`npar',0)
	loc pp = 3-`gfix'
	loc i 1
	while `i'<=_N {
		loc uze = `touse'[`i']
		if `uze' {
* multiply centile derivs by data for each parameter
			loc param 1
			loc k 0
			while `param'<=`pp' {
				scalar `h' = `h`param''[`i']
				loc np `n`param''
				loc j 0
				while `j'<`np' {
					loc j = `j'+1
					loc k = `k'+1
					if `j'<`np' {
						scalar `xx'=`xn`k''[`i']
					}
					else scalar `xx'=1
					matrix `hbeta'[1,`k'] =`h'*`xx'
				}	
				loc param = `param'+1
			}
			matrix `hVh' = `hbeta'*`V'
			matrix `hVh' = `hVh'*`hbeta''
			scalar `se'=sqrt(`hVh'[1,1])
			if `no_x' { loc i=_N }
			else qui replace `cse'=`se' in `i'
		}
		loc i = `i'+1
	}
	if `no_x' { qui replace `cse' = `se' if `touse' }
end

program define _se4 /* SE for estimated centiles, MPN and MEN models */
	loc dnum  `1'	/* distribution number */
	loc touse `2'	/* filter for original obs */
	loc vars1 `3'	/* covariates for M-curve */
	loc vars2 `4'	/* covariates for S-curve */
	loc vars3 `5'	/* covariates for G-curve */
	loc vars4 `6'	/* covariates for D-curve */
	loc base1 `7'	/* base vars  for M-curve */
	loc base2 `8'	/* base vars  for S-curve */
	loc base3 `9'	/* base vars  for G-curve */
	loc base4 `10'/* base vars  for D-curve */
	loc M `11'
	loc S `12'
	loc G `13'
	loc D `14'
	loc V `15'	/* VCE matrix of parameters */
	loc cv = "`16'"!=""
	loc tau = "`17'"!=""
	loc lns = "`18'"!=""
	loc gfix = "`19'"!=""
	loc dfix = "`20'"!=""
	loc centile `21'
	loc cse `22'	/* to hold SE, assumed to exist */
	loc no_x `23'
* Put names of covariates into xn*

	loc small 1e-6

	loc p 1	/* parameter no.: 1=m, 2=s, 3=g, 4=d. */
	loc npar 0	/* npar counts the covariates, including _cons's */
	while `p'<=4 {
		loc no_gpar = `gfix' & (`p'==3)
		loc no_dpar = `dfix' & (`p'==4)
		if `no_gpar' | `no_dpar' { loc n`p' 0 }
		else loc n`p' 1
*noi di in red "p=" `p' ", np=" `n`p'' ", varsp= |`vars`p''|, basep= |`base`p''|"
		if "`vars`p''`base`p''"!="" {
/*
	Put names of the FP and base covariates into xn*
*/
			parse "`vars`p'' `base`p''", parse(" ")
			while "`1'"!="" {
				loc n`p'=`n`p''+1
				loc npar=`npar'+1
				loc xn`npar' "`1'"
				mac shift
			}
		}
/*
	Add one param for _cons
*/
		if !`no_gpar' & !`no_dpar' { loc npar=`npar'+1 }
		loc p=`p'+1
	}
/*
	Calc first derivatives w.r.t. params.
*/
	tempvar h1 h2 h3 h4 A B DD zq s
	cap confirm var `centile'
	if _rc {
		tempname uq
		scalar `uq'=invnorm(`centile'/100)
	}
	else local uq `centile' /* uq contains Z-scores (zci option) */
	gen `DD'=1+`D'*abs(`uq')
	gen `zq'=sign(`uq')*(`DD'^(1/`D')-1)
	if `lns' {	/* log S parameterization */
		gen `s'=exp(`S')
	}
	else gen `s'=`S'
	if `cv' { replace `s'=`s'*`M' } /* `s' is sigma, not cv */

	if `dnum'==7 {	/* MPN */
	    if `tau' {		/* tau parameterization */
		tempvar t	/* lambda */
		gen `t'=1-`G'*`M'/`s'
	    }
	    else loc t `G'
	    gen `B'=`zq'*`s'/`M'
	    gen `A'=cond(abs(`t')<`small', exp(`B'), (1+`t'*`B')^(1/`t'))
	    gen `h1'=cond(abs(`t')<`small', `A'*(1-`B'), `A'*(1-`B'/`A'^`t'))
	    if `cv' { gen `h2'=`M'*`zq'*`A'^(1-`t') }
	    else gen `h2'=`zq'*`A'^(1-`t')
	    if !`gfix' {
		gen `h3'=cond(abs(`t')<`small', -0.5*`M'*`A'*`B'^2, /*
		 */ (`M'*`A'/`t')*(`B'/`A'^`t'-log(`A')) )
	    }
	    if !`dfix' {
		gen `h4'=(`zq'+sign(`uq'))*(`s'/`D')*`A'^(1-`t')* /*
		 */ (abs(`uq')/`DD'-log(`DD')/`D')
	    }
	}
	else if `dnum'==8 { /* MEN */
	    gen `A'=1+`G'*`zq'
	    gen `h1'=1
	    if `cv' {
	        gen `h2'=cond(abs(`G')<`small', `M'*`zq', /*
		*/ `M'*log(`A')/`G' )
	    }
	    else gen `h2'=cond(abs(`G')<`small', `zq', log(`A')/`G' )
	    if !`gfix' {
		gen `h3'=cond(abs(`G')<`small', /*
		 */ -0.5*`s'*`zq'^2, (1-1/`A'-log(`A'))*`s'/`G'^2 )
	    }
	    if !`dfix' {
		gen `h4'=(`zq'+sign(`uq'))*(`s'/(`D'*`A'))* /*
		 */ (abs(`uq')/`DD'-log(`DD')/`D')
	    }
	}
* Calc SE of centile, element-by-element
	tempname xx h hbeta hVh se
	matrix `hbeta'=J(1,`npar',0)
	loc pp=4-`gfix'-`dfix'
	loc i 1
	while `i'<=_N {
		loc uze=`touse'[`i']
		if `uze' {
* multiply centile derivs by data for each parameter
			loc param 1
			loc k 0
			while `param'<=`pp' {
				scalar `h'=`h`param''[`i']
				loc np `n`param''
*noi di in red "obs i=`i', param=`param', np=`np', h=" `h'
				loc j 0
				while `j'<`np' {
					loc j=`j'+1
					loc k=`k'+1
					if `j'<`np' {
						scalar `xx'=`xn`k''[`i']
					}
					else scalar `xx'=1
*noi di in red "j=`j', k=`k', xx=" `xx'
					matrix `hbeta'[1,`k']=`h'*`xx'
				}	
				loc param=`param'+1
			}
			matrix `hVh'=`hbeta'*`V'
			matrix `hVh'=`hVh'*`hbeta''
			scalar `se'=sqrt(`hVh'[1,1])
			if `no_x' { loc i=_N }
			else qui replace `cse'=`se' in `i'
		}
		loc i = `i'+1
	}
	if `no_x' { qui replace `cse' = `se' if `touse' }
end

*! version 1.0.1  07/17/93
* Stata factory program mx_out, modified to remove model chisquare output.
program define mx_out
	parse "`*'", parse(" ,")
	local cmd `1'
	mac shift
	local options "*"
	parse "`*'"
	if ("`cmd'"!="$S_E_cmd") { error 301 }
	di
	di in gr "$S_E_ttl" _col(53) "Number of obs    =" in yel %8.0f $S_E_nobs
	di in gre "Log Likelihood =" in yel %15.7f $S_E_ll _c
	if ("$S_E_pr2"!="") { 
		di _col(22) in gr "Pseudo R2        =" /*
			*/ in yel %8.4f $S_E_pr2 _c 
	}
	di _n

	matrix mlout, `options'
end