;+
; NAME:
; get_w
;
;
; PURPOSE:
; compute the "w" parameter for boxcox (utility routine)
;
;
; INPUTS:
; lambda: lambda parameter from boxcox
;
;
; OUTPUTS:
; the w parameter from boxcox
;
;
;-
FUNCTION get_w, lambda
  COMMON boxcox, x1, y1
  k2 = geo_mean(y1, /double)
  
  IF lambda NE 0 THEN BEGIN
     k1 = 1./(lambda*k2^(lambda - 1.))
     w = k1 * (y1^lambda - 1.)
  ENDIF ELSE BEGIN
     w = k2 *alog(y1)
  ENDELSE 
  
  sse = sse(w)
  RETURN, sse
END

;+
; NAME:
; boxcox
;
;
; PURPOSE:
; determine the Cox transformation that makes the data most normal
;
; from Applied linear statistical models
; Neter, Kutner, et al.
; pp 133
;
; CATEGORY:
; Stats
;
;
; INPUTS:
; x: the 'x' values
; y: the 'y' values (must be the same length as x)
;
;
; KEYWORD PARAMETERS:
; DOUBLE: do the calculation in double precision
; PLOT: plot a diagnostic plot
; RANGE: the range of Cox transformations to search for the answer,
; default [-2,2]
; SUBTITLE: set the subtitle for the diagnostic plot
;
;
; OUTPUTS:
; structure containing:
;  'minsse' - the minimum sse
;  'lambda' - the Cox transformation where that is true
;
; COMMON BLOCKS:
; BOXCOX: contains the data passed in to then make it availiable to
; the utility routine get_w
;
;
;
; RESTRICTIONS:
; all the Y values must be positive and non-zero
;
;
; EXAMPLE:
; x=findgen(10)    
; y=randomn(seed, 10)
; y-=min(y)-0.1
; ans = boxcox(x, y) 
; help, ans, /str
; ** Structure <1eaa814>, 2 tags, length=8, data length=8, refs=1:
;    MINSSE          FLOAT           4.78825
;    LAMBDA          FLOAT     Array[1]
; print, ans.lambda
;     0.700000
;
; MODIFICATION HISTORY:
;
;       Tue Jan 22 14:11:45 2008, Brian Larsen
;		documented and fixed up
;
;-
FUNCTION boxcox, x, y, double=double, PLOT=plot, RANGE=range, SUBTITLE=subtitle
  COMMON boxcox, x1, y1
  x1=x
  y1=y
  
  IF N_ELEMENTS(range) NE 2 THEN range=[-2,2]
  
  ans = fltarr(41)
  lams = (findgen(41)-20)/10.
  FOR i=0l,N_ELEMENTS(ans)-1  DO BEGIN
     ans[i] = get_w( (i-((N_ELEMENTS(ans)-1)/2.) )/10. )
  ENDFOR
  
  ind = where(lams LT range[0] OR lams GT range[1], complement=comp)
  lams = lams[comp]
  ans  = ans[comp]
  
  IF KEYWORD_SET(plot)  THEN BEGIN
     plot, lams, ans, psym=-1, $
           title='BOXCOX', $
           xtitle='BOXCOX power trans', $
           ytitle='SSE', $
           subtitle=subtitle
     plots, lams[where(ans EQ min(ans))], ans[where(ans EQ min(ans))], psym=2, symsize=3
  ENDIF
  
  ans = create_struct('minsse', min(ans), $
                      'lambda', lams[where(ans EQ min(ans))] )
  
  RETURN, ans
END