%macro manly(phimin=,phimax=, steps=, model=, class=, stmts=, data=, response=);

/***********************************************************/
/* This macro performs the exponential transformation by   */
/* Manly, B. F. J. (1976) Exponential data transformations.*/
/* Statistician, 25, 37-42.                                */
/*                                                         */
/*   phimin   = smallest phi                               */
/*   phimax   = largest phi                                */
/*   steps    = number of steps in grid search             */
/*   class    = list of variables for class statment       */
/*   model    = list of terms in model statement           */
/*   stmts    = %str(other mixed statements, incl. random) */
/*   data     = dataset with original data                 */
/*   response = name of response variable in dataset       */
/*                                                         */
/***********************************************************/
data t;
set &data;
phi=&phimin;
size=(&phimax-&phimin)/&steps;
do i=0 to &steps;
  if abs(phi)>1e-10 then _y_=(exp(phi*&response)-1)/phi;
  if abs(phi)<1e-10 then _y_=&response;
  jac=phi*&response;
  output;
  phi=phi+size;
end;

proc sort data=t out=t;
by phi;
run;

*proc print data=t; run;

ods output fitstatistics=fit;
proc mixed data=t method=ml;
class &class;
model _y_=&model;
&stmts;
by phi;
run;

data fit2;
set fit;
if descr='-2 Log Likelihood';
loglik=-value/2;

proc means data=t;
var jac;
output out=jac sum=;
by phi;
run;

data loglik;
merge jac fit2;
loglik=loglik+jac;

proc print data=loglik;
run;

symbol value=dot i=join;
proc gplot data=loglik;
plot loglik*phi;
run;

proc means data=loglik;
var loglik;
id phi;
output out=max max=max;
run;

data r;
if _n_=1 then set max;
set loglik;
if loglik=max;

title 'optimal value for phi';
proc print data=r;
var phi loglik;
run;
title ' ';

%mend;


data s;
do block=1 to 50;
  b=ranexp(12);
  do trt=1 to 4;
     x=b  + ranexp(1);
	 output;
  end;
end;
;
proc print data=s; run;

proc univariate;
histogram x;
run;

%manly(phimin=-1,phimax=1, steps=10, model=, class=trt block, data=s, stmts=%str(random block;parms (1)(1);), response=x);