A Robust Measure of Scale modeled on Maronnaa & Zamarb's performance improvements on Rousseeuw and Croux's skewness and efficiency improvements on median absolute deviation. Blatantly transposed from the R robustbase module.

robustTau.py

from scipy import stats as scistats
import numpy as np

# Implementation of Tau from http://amstat.tandfonline.com/doi/abs/10.1198/004017002188618509#.VDgKhdR4rEh
# blatently transposed R robustbase library from http://r-forge.r-project.org/scm/?group_id=59, OGK.R
def scaleTau2(x, c1 = 4.5, c2 = 3.0, consistency = True, mu_too = False, *xargs, **kargs):
	## NOTA BENE: This is *NOT* consistency corrected
	x = np.asarray(x)
	n = len(x)
	medx = np.median(x) 
	x_abs = abs(x - medx)
	sigma0 = median(x_abs)
	if c1 > 0:
		## w <- pmax(o, 1-(x_abs / (sigma0 * c1))^2)^2 -- but faster: 
		x_abs = x_abs / (sigma0 * c1)
		w = 1 - x_abs * x_abs
		w = ((abs(w) + w) /2)**2
		mu = sum(x * w) / sum(w) # Gabe: implicit float?
	else:
		mu = medx

	x = (x - mu) / sigma0
	rho = x**2
	rho[rho > c2**2] = c2**2
	## sigma2 <- sigma0**2 * sum(rho) / 2

	if consistency:
		def Erho(b):
			## E [ rho_b ( X ) ]    X ~ N(0,1)
			# GABE: norm.cdf = pnorm, norm.pdf = dnorm.  based on http://adorio-research.org/wordpress/?p=284
			return 2*((1-b**2) * scistats.norm.cdf(b) - b * scistats.norm.pdf(b) + b**2) - 1
		def Es2(c2):
			## k^2 * E[ rho_{c2} (X' / k) ] , where X' ~ N(0,1), k = qnorm(3/4)
			return Erho(c2 * scistats.norm.ppf(3/float(4)))
		## the asymptotic E[ sigma**2(X) ] is Es2(c2):
		## TODO: 'n-2' below will probably change; therefore not yet documented
		if consistency == "finiteSample":
			nEs2 = (n-2) * Es2(c2)
		else:
			nEs2 = n * Es2(c2)
	else:
		nEs2 = n

	## sqrt(sigma2) == sqrt( sigma0^2 / n * sum(rho) ) :
	tau = sigma0 * sqrt(sum(rho)/nEs2)

	if mu_too:
		return mu, tau
	else:
		return tau