*! v 3.1.0 11-Sep-94.
program define _jprfpfp
version 3.1
loc anova `1'
loc blogit `2'
loc cmd `3'
loc continu "`4'"
loc devthr "`5'"
loc df `6'
loc df3 `7'
loc df4 `8'
loc dist `9'
loc exp "`10'"
loc fast "`11'"
loc fixpowe "`12'"
loc glm `13'
loc gplot `14'
loc graph `15'
loc if `16'
loc in `17'
loc kx "`18'"
loc lin `19'
loc log "`20'"
loc lnx `21'
loc meq2 `22'
loc mnlnwt `23'
loc name "`24'"
loc nobs `25'
loc normal `26'
loc np `27'
loc options "`28'"
loc origin "`29'"
loc p `30'
loc pwrs "`31'"
loc qreg `32'
loc repeat "`33'"
loc rhs "`34'"
loc saving "`35'"
loc search `36'
loc shift "`37'"
loc x `38'
loc y `39'
loc zero "`40'"
loc zeta "`41'"
loc small 1e-6
if `normal' {
	loc devcal "\$_nobs*(1+log(2*_pi*_result(4)/\$_nobs)-\$_mnlnwt)"
	loc pcal "fprob(\$_n1,\$_n2,(exp(\$_d/\$_nobs)-1)*(\$_n2)/(\$_n1))"
}
else {
	loc pcal "chiprob(\$_n1,\$_d)"
	if `dist'==2 { loc devcal "-2*\$S_E_ll" }
	else if `glm' {
		gl S_3 1
		loc devcal "\$S_3"
	}
	else loc devcal "-2*_result(2)"
}
if `search' {
	cap est drop fpres
	cap est drop temp
	cap conf num $FP_nmod
	if !_rc {
		loc j 0
		while `j'<$FP_nmod {
			loc j=`j'+1
			gl FP_pw`j'
			gl FP_dv`j'
			gl FP_ga`j'
		}
		gl FP_nmod
	}
	gl FP_fixp
	loc nm 0
	tempvar xp1 xp2 xp3
	quietly {
		gen `xp1'=.
		gen `xp2'=.
		if "`devthr'"!="" {
			tempvar yl yh stdp
			gen `yl'=.
			gen `yh'=.
			gen `stdp'=.
		}
	}
	loc base "$S_E_base"
	if "`fixpowe'"!="" {
		tempvar xpp
		qui gen `xpp'=.
		_jprfppp `xpp' "`fixpowe'" 0 0
		if $S_2 {
			di in red "fixpowers must be in increasing order"
			exit 198
		}
		loc nxp $S_1
		loc j 1
		loc h
		while `j'<=`nxp' {
			loc hj "h`j'"
			tempvar `hj'
			qui gen ``hj''=.
			loc h "`h' ``hj''"
			loc j=`j'+1
		}
		_jprfpgn `x' `lnx' `xpp' `nxp' `h'
/*
	Determine positions of terms in `fixpowe'
*/
		loc j 1
		while `j'<=`np' {
			loc pj=`p'[`j']
			loc i 1
			loc pos 0
			while `i'<=`nxp' {
				if abs(`xpp'[`i']-`pj')<`small' {
					loc pos `i'
				}
				loc i=`i'+1
			}
			loc pos`j' `pos'
			loc j=`j'+1
		}
	}
	else {
		loc nxp 0
		loc j 1
		while `j'<=`np' {
			loc pos`j' 0
			loc j=`j'+1
		}
	}
/*
	Calc deviance=-2(log likelihood) for base regression,
	allowing for possible weights.
*/
	if `anova' {
		if "`continu'"!="" { loc cts "cont(`continu')" }
	}
	qui `cmd' `y' `base' `exp' if `x'!=., `cts' `options'
	loc rdf=_result(5)
	if `rdf'<(1+`df'+`nxp') { error 2001 }
	if `qreg' { loc dev0=`nobs'*log($S_E_msd/$S_E_rsd) }
	else loc dev0=`devcal'
	loc devm2 `dev0'
	loc devm1 `dev0'
	if "`log'"!="" {
		di in gr "Power1" _col(10) "Power2" /*
		 */ _col(21) "Dev diff" _col(35) _cont
		if `normal' { di in gr "Res SD" }
		else di _n
	}
	qui if `graph' {
		loc maxd=`np'+`meq2'*(`np'*(`np'+1))/2
		if `maxd'>_N { set obs `maxd' }
		tempvar p1 p2 dv
		gen `p1'=.
		gen `p2'=.
		gen `dv'=.
		loc g 1
	}
	loc j 1
	while `j'<=`np' {
		loc pj=`p'[`j']
		loc pk .
		if `pos`j''>0 {
			qui replace `xp1'=`h`pos`j'''*`lnx'
		}
		else {
			if abs(`pj')<`small' {
				qui replace `xp1'=`lnx'
			}
			else {
				qui replace `xp1'=cond(`x'==0,0,`x'^`pj')
			}
		}
		if !`df4' {
			if `anova' { loc cts "cont(`continu' `xp1' `h')" }
			qui `cmd' `y' `xp1' `h' `base' `exp', `cts' `options'
			if `qreg' { loc dev=`nobs'*log($S_E_msd/$S_E_rsd) }
			else loc dev=`devcal'
			loc nm=`nm'+1
			gl FP_pw`nm' `pj'
			gl FP_dv`nm' `dev'
			gl FP_ga`nm'=_b[`xp1']
			if `normal' { loc rsd=_result(9) }
			_jprfpdt "`devthr'" `dev' `pj' `pk' `yl' `yh' `stdp'
			qui if `graph' & `gplot'<2 {
				replace `p1'=`pj'+.05 in `g'
				replace `dv'=`dev' in `g'
				loc g=`g'+1
			}
			if `dev'<`devm1' {
				loc devm1 `dev'
				loc pm1 `pj'
				loc rsd1 `rsd'
				if !`meq2' { gl FP_best `nm' }
			}
			if "`log'"!="" {
				di in ye %6.3f `pj' _col(12) "." _col(20) /*
				 */  %9.3f `dev0'-`dev' _col(32) _cont
				if `normal' { di " " %8.0g `rsd' }
				else di
			}
		}
		if `meq2' {
			if "$FP_BOTH"=="" { loc k = `j'+("`repeat'"!="") }
			else loc k 1
			loc k1 `k'
			loc devdif 0
			while `k'<=`np' {
			 if `df3' {
				loc k `np'
				loc pk 1
			 }
			 else loc pk=`p'[`k']
			 if `pk'==`pj' { /* Case p1=p2 */
			 	qui replace `xp2'=`xp1'*`lnx'
			 }
			 else {
			 	if `pos`k''>0 {
			 		qui replace `xp2'=/*
					 */ `h`pos`k'''*`lnx'
			 	}
			 	else {
			 		if abs(`pk')<`small' {
			 			qui replace `xp2'=`lnx'
			 		}
			 		else {
	 			qui replace `xp2'=cond(`x'==0,0,`x'^`pk')
			 		}
			 	}
			 }
			 if `anova' {
			 	loc cts "cont(`continu' `xp1' `xp2' `h')"
			 }
			 qui `cmd' `y' `xp1' `xp2' `h' /*
			  */ `base' `exp', `cts' `options'
			 if `qreg' { loc dev=`nobs'*log($S_E_msd/$S_E_rsd) }
			 else loc dev=`devcal'
			 if `normal' { loc rsd=_result(9) }
			 loc nm=`nm'+1
			 gl FP_pw`nm' "`pj',`pk'"
			 gl FP_dv`nm' `dev'
			 gl FP_ga`nm'=_b[`xp2']/_b[`xp1']
			 _jprfpdt "`devthr'" `dev' `pj' `pk' `yl' `yh' `stdp'
			 qui if `graph' & `gplot'!=1 {
				replace `p1'=`pj' in `g'
				replace `p2'=`pk' in `g'
				replace `dv'=`dev' in `g'
				loc g=`g'+1
			 }
/*
	Increasing deviance is devdif<0.
*/
			 if `k'>`k1' { loc devdif=`devold'-`dev' }
			 loc devold `dev'
			 if `dev'<`devm2' {
			 	loc devm2 `dev'
			 	loc pm21 `pj'
			 	loc pm22 `pk'
			 	loc rsd2 `rsd'
				gl FP_best `nm'
			 }
			 if `df4'&`dev'<`devm1'&(`pj'==1|`pk'==1) {
				loc devm1 `dev'
				if `pj'==1 { loc pm1 `pk' }
				else loc pm1 `pj'
				loc rsd1 `rsd'
			 }
			 if "`log'"!="" {
			 	di in ye %6.3f `pj' _col(10) /*
			 	 */ %6.3f `pk' _col(20) /*
			 	 */ %9.3f `dev0'-`dev' _col(32) _cont
			 	if `normal' {
			 		di in ye " " %8.0g `rsd'
			 	}
			 	else di
			 }
/*
	If deviance is increasing, abort sequence
*/
			 if "`fast'"=="fast" & `devdif'<0 {
			 	loc k `np'
			 }
			 loc k=`k'+1
			}
		}
		loc j=`j'+1
	}
	if `anova' { loc cts "cont(`continu' `x')" }
	qui `cmd' `y' `x' `base' `exp', `cts' `options'
	if `qreg' { loc devlin=`nobs'*log($S_E_msd/$S_E_rsd) }
	else loc devlin=`devcal'
	qui if !`lin' {
		replace `xp2'=`x'^2
		if `anova' { loc cts "cont(`continu' `x' `xp2')" }
		`cmd' `y' `x' `xp2' `base' `exp', `cts' `options'
		if `qreg' { loc devquad=`nobs'*log($S_E_msd/$S_E_rsd) }
		else loc devquad=`devcal'
		gen `xp3'=`x'*`xp2'
		if `anova' { loc cts "cont(`continu' `x' `xp2' `xp3')" }
		cap `cmd' `y' `x' `xp2' `xp3' `base' `exp', `cts' `options'
		if `qreg' { loc devcub=`nobs'*log($S_E_msd/$S_E_rsd) }
		else loc devcub=`devcal'
		drop `xp3'
/*
	BoxTid model: x and x*ln x
*/
		replace `xp2'=`x'*`lnx'
		if `anova' { loc cts "cont(`continu' `x' `xp2')" }
		`cmd' `y' `x' `xp2' `base' `exp', `cts' `options'
		if `qreg' { loc devBT=`nobs'*log($S_E_msd/$S_E_rsd) }
		else loc devBT=`devcal'
	}
/*
	Tests (F or chi-square) based on deviances.
*/
	#delimit ;
	loc d0lin=`dev0'-`devlin'; loc d `d0lin';
	loc n1 1; loc n2=`rdf'-1;
	loc P0lin=`pcal';

	loc dlinm1=`devlin'-`devm1'; loc d `dlinm1';
	loc n1=1+`nxp'+`df4'; loc n2=`rdf'-1-`n1';
	loc Plinm1=`pcal';

	loc d0m1=`dev0'-`devm1'; loc d `d0m1';
	loc n1=2-`lin'+`nxp'+`df4'; loc n2=`rdf'-`n1';
	loc P0m1=`pcal';

	if !`lin' {
		loc d1quad=`devlin'-`devquad'; loc d `d1quad';
		loc n1 1; loc n2=`rdf'-2;
		loc P1quad=`pcal';

		loc dlinBT=`devlin'-`devBT'; loc d `dlinBT';
		loc PlinBT=`pcal';

		loc d2cub=`devquad'-`devcub'; loc d `d2cub';
		loc n2=`rdf'-3;
		loc P2cub=`pcal';
	};
	if `meq2' {
		loc rdf4=`rdf'-4-`nxp'+`df3';
		loc dm1m2=`devm1'-`devm2'; loc d `dm1m2';
		loc n1=2-`df3'-`df4'; loc n2 `rdf4';
		loc Pm1m2=`pcal';
		loc d0m2=`dev0'-`devm2'; loc d `d0m2';
		loc n1=4+`nxp'-`df3';
		loc P0m2 =`pcal';
		loc dlinm2=`devlin'-`devm2'; loc d `dlinm2';
		loc n1=3+`nxp'-`df3';
		loc Plinm2 =`pcal';
		loc dquadm2=`devquad'-`devm2'; loc d `dquadm2';
		loc n1=2+`nxp'-`df3';
		loc Pquadm2 =`pcal';
	};
	#delimit cr
	if !`lin' {
		if `graph' {
			qui replace `dv'=`devlin'-`dv'
			graph `dv' `p1', xla(`pwrs') yla s([`p2']) /*
			*/ l1(Gain) b2("Power") xli(`pwrs') `saving'
		}
		drop `p' `lnx' `xp1' `xp2'
	}
	if `meq2' { loc pwrs "`pm21' `pm22' `fixpowe'" }
	else loc pwrs "`pm1' `fixpowe'"
	if "`devthr'"!="" {
		cap drop _etal
		cap drop _etah
		cap drop _stdp
		rename `yl' _etal
		rename `yh' _etah
		rename `stdp' _stdp
	}
}
/*
	Re-estimate final (p~) model if search occurred,
	otherwise estimate fixed-powers model.
*/
if "`name'"=="" { loc name X }
cap drop `name'
rename `x' `name'
lab var `name' "`rhs'"
fpgen `name', powers(`pwrs') replace `zero'
loc xp "$S_1"
if `anova' { loc cts "cont(`continu' `xp')" }
qui `cmd' `y' `xp' `base' `exp', `cts' `options'
if !`search' | `lin' {
	if `qreg' { loc dev0=`nobs'*log($S_E_msd/$S_E_rsd) }
	else loc dev0=`devcal'
}
else {
	gl FP_curr $FP_best
	gl FP_nmod `nm'
	gl FP_fixp "`fixpowe'"
}
/*
	Next line is a kludge since can't pass
	long string (`base') as argument.
*/
gl S_E_base "`base'"
#delimit ;
_jprfpse `blogit' `cmd' "`continu'" "`d0lin'" "`d0m1'" "`d0m2'" "`d1quad'"
 "`d2cub'" `dev0' "`devBT'" "`devcub'" "`devlin'" "`devm1'" "`devm2'"
 "`devquad'" `df' `dist' "`dlinBT'" "`dlinm1'" "`dlinm2'" "`dm1m2'"
 "`dquadm2'" "`exp'" "`fixpowe'" "`if'" "`in'" "`kx'" `lin' `meq2' "`name'"
 `nobs' "`nxp'" "`options'" "`origin'" "`P0lin'" "`P0m1'" "`P0m2'" "`P1quad'"
 "`P2cub'" "`PlinBT'" "`Plinm1'" "`Plinm2'" "`pm1'" "`Pm1m2'" "`pm21'"
 "`pm22'" "`Pquadm2'" "`pwrs'" "`rhs'" "`rsd1'" "`rsd2'" `search' "`shift'"
 "`xp'" `y' "`zero'" "`zeta'" ;
#delimit cr
_jprfprp, summary `estimat' `compari'
end