*
*  BOXCOX.PRG
*  Examples from Sections 7.10 and 7.11
*
all 200
*
*   Create data with a square root transformation.
*
set x = %rangamma(5.0)
set sqry = .5 + 3*sqrt(x) + %ran(4.0)
set y = sqry**2
*
*   Using MAXIMIZE - from section 7.10. The free coefficients
*   b0, b1 and sigmasq will be transforms of the .5, 3 and 16=4**2
*   used in generating the data since sqrt(y) and sqrt(x) are
*   translated, scaled versions of the Box-Cox transformed data.
*
nonlin b0 b1 lambda sigmasq
frml rhsfrml  = b0+b1*%boxcox(x,lambda)
frml boxcox   = (lambda-1)*log(y) + %logdensity(sigmasq,%boxcox(y,lambda)-rhsfrml(t))
*
*   This estimates a simple linear regression (special case when l=1),
*   and sets the initial conditions based upon it.
*
linreg y
# constant x
compute sigmasq=%seesq,b0=%beta(1)+%beta(2)-1
compute b1=%beta(2),lambda=1.0
maximize(method=bfgs) boxcox
*
*   Using FIND - from section 7.11
*
compute sumlogy=%sum(%log(y))
nonlin lambda
compute lambda=1.0
declare real loglikely
find maximum loglikely
   set xbc = %boxcox(x,lambda)
   set ybc = %boxcox(y,lambda)
   linreg(noprint) ybc
   # constant xbc
   compute sigmasq=%rss/%nobs
   compute loglikely=(lambda-1)*sumlogy-.5*%nobs*log(sigmasq)
end find