;+
; 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