********************************************************************************************************************************
* pbreg fits the Preece and Baines 1978 family of growth curves.
* version 1.0
* Author: Adrian Sayers
* Date: 06.03.2013
*
********************************************************************************************************************************

prog define pbreg , rclass
version 9.2
syntax varlist(min=2 max=2 numeric) [if] [in] , 	[ 		Model(integer 1)		///
															h1( real 0)				///
															HTheta(real 0)			///
															s0(real 0)				///
															s1(real 0)				///
															THeta(real 0)			///
															GAmma(real 1)			///
															p0(real 0)				///
															p1(real 0)				///
															q1(real 0)				///
															level(real 95) 			///
															nolog					///
															trace					///
															iterate(integer 10000)  ///
															grid					///
															GRIDSPace(integer 5)    ///
															bestvals(integer 5) 	///
															cdif					///
															eps(real 1e-5)			///
													]

tokenize `varlist'
	tempname outcome time
		gen `outcome' = `1'
		gen `time' = `2'

*********************************************************************************************
* Speed Model convergence, zero data, starting values and grid search
*********************************************************************************************
quietly { // find initial values.

if "`cdif'"!="" { // Calculate first central difference derivative as starting values.
	tempvar vel

	sort `time'
		gen `vel' = (`outcome'[_n+1] - `outcome'[_n-1]) / (`time'[_n+1] -`time'[_n-1])
			su `vel' , meanonly
				local phv =r(max)
			su `outcome' , meanonly
				local h1  = r(max)
				su `outcome' if `vel'>`phv'-0.000001  & `vel'<`phv'+0.000001 , meanonly
					local htheta = r(mean)
				su `time' if `vel'>`phv'-0.000001  & `vel'<`phv'+0.000001 , meanonly
					local theta = r(mean)
					local s0 = 1
					local s1 = 0.1
				} // end cdif starting values.


if "`grid'" != "" { // Grid Search for starting values.
	tempfile start_values

	cap postclose sv
		postfile sv rss h1 htheta theta using `start_values' , replace

	// Use tokens
	 su `outcome'
		local min_ht = r(min)
		local max_ht = r(max)
		local grid_ht = (r(max) - r(min))/`gridspace'
	 su `time'
		local min_age = r(min)
		local max_age = r(max)
		local grid_age  = (r(max) - r(min))/`gridspace'
	noisily di "Conducting grid search"
	forvalues h1 = `min_ht'(`grid_ht')`max_ht'  {
		forvalues htheta = `min_ht'(`grid_ht')`max_ht'  {
			forvalues theta = `min_age'(`grid_age')`max_age'  {
				capture pbreg `outcome' `time' , model(1) h1(`h1') htheta(`htheta') theta(`theta') s0(1) s1(0) iterate(1)
					post sv (`e(rss)' ) ( `h1') (`htheta') (`theta')
					 noisily di "*" _c
			}
		}
	}
	postclose sv

	// Store the best 5 starting values
		preserve
			use  `start_values' , clear
				sort  rss
				 noisily di _n "Preserving best starting values" _n
					forvalues  i = 1 / `bestvals' {
						local h1_`i'		= h1[`i']
							local htheta_`i' 	= htheta[`i']
								local theta_`i'		= theta[`i']
									noisily di "*" _c
										}
		restore

	// Try 5 of the best to completeion
	local minssr = .
	local best_h1 = .
	local best_htheta = .
	local best_theta = .


noisily di _n "Trying best starting values" _n
	forvalues i = 1 / `bestvals' {
	capture  pbreg `outcome' `time' , model(1) h1(`h1_`i'') htheta(`htheta_`i'') theta(`theta_`i'') s0(1) s1(0) iterate(`iterate')
		if e(rss)  <`minssr' {
			local best_h1 = `h1_`i''
			local best_htheta = `htheta_`i''
			local best_theta = `theta_`i''
			local minssr = e(rss)
							} // end if rss is smaller
					noisily di "*" _c
						}

	// Initial Values

	local h1		=_b[/h1]
	local htheta    =_b[/htheta]
	local theta		=_b[/theta]
	local s0        =_b[/s0]
	local s1        =_b[/s1]
 } // End Grid Search

 } // end quietly find iinitial values.

if `model'==1 local initial "h1 `h1'  htheta `htheta' s0 `s0' s1 `s1' theta `theta'"
if `model'==2 local initial "h1 `h1'  htheta `htheta' s0 `s0' s1 `s1' theta `theta' ga `gamma'"
if `model'==3 local initial "h1 `h1'  htheta `htheta' P0 `p0' P1 `p1' theta `theta' Q1 `q1'"


*********************************************************************************************
* Fit the Models
*********************************************************************************************
capture noisily  {
*********************************************************************************************
// Model 1
*********************************************************************************************

if `model' ==1 {
	 nl (`outcome' = {h1} -  ((2*({h1}-{htheta})) / ( exp({s0}*(`time'-{theta})) + exp({s1}*(`time'-{theta}))))) `if' `in' , ///
		initial(`initial' ) level(`level') iterate(`iterate') eps(`eps')   `log' `trace'
				} // end model 1

*********************************************************************************************
// Model 2
*********************************************************************************************

if `model' ==2 {
nl 	(`outcome' = {h1} - (({h1}-{htheta}) / (((0.5*exp(({ga}*{s0})	*(`time'-{theta}))) + (0.5*exp(	({ga}*{s1}) * 	(`time'-{theta}))))^(1/{ga})))) `if' `in' , initial(`initial' )  level(`level')  iterate(`iterate')
				} // end model 2

*********************************************************************************************
// Model 3
*********************************************************************************************

if `model' ==3 {
nl 	(`outcome' = {h1} - ((4*({h1}-{htheta})) /  ((exp({P0}*(`time'-{theta})) +  exp({P1}*(`time'-{theta})))* (1+ exp({Q1}*(`time'-{theta})))))) `if' `in' , initial(`initial' ) level(`level')
       			} // end model 3

       			}
if _rc != 0 {
	 di as err "PB Model Failure "
	 exit _rc
	 		}


end