moments.jl

# * Fri Feb 24 2012 09:09:10 PM Steven E. Pav <shabbychef@gmail.com>
#
# computation of moments and moment-based statistics
# 
# see also: 
# * http://www.johndcook.com/blog/2008/09/26/comparing-three-methods-of-computing-standard-deviation/
# * http://www.johndcook.com/standard_deviation.html
# * J. Bennett, et. al., 'Numerically Stable, Single-Pass,
#   Parallel Statistics Algorithms,' Proceedings of IEEE
#   International Conference on Cluster Computing, 2009.
#   http://www.janinebennett.org/index_files/ParallelStatisticsAlgorithms.pdf
# * T. Terriberry, 'Computing Higher-Order Moments Online,' 
#   http://people.xiph.org/~tterribe/notes/homs.html

# compute the kth moment of a vector
# based on the divide and conquer scheme from
# Bennett et. al.

# 2FIX: make this private?
# and not slow as hellll.
function rec_moments(v::AbstractVector, ord::Int)
  # returns the number of elements, the mean, and 
	# a (ord - 1)-vector consisting of the
	# 2nd through ord'th centered sum, defined
	# as sum_j (v[j] - mean)^i
	# for i > 1
	n = length(v)
	if n == 1
		return (n, v[1], zeros(eltype(v),ord-1,1))
	elseif n == 2  # a premature optimization not really needed.
		mu = 0.5 * sum(v)
		moms = zeros(eltype(v),ord-1,1)
		for pp=2:ord
			moms[pp-1] = (v[1] - mu) ^ (convert(Float,pp)) + (v[2] - mu) ^ (convert(Float,pp))
		end
		return (n, mu, moms)
	else
		rets = zeros(eltype(v),ord,1)
		cut = convert(Int,floor(0.5 * n))
		(n1,mu1,moms1) = rec_moments(v[1:cut],ord)
		(n2,mu2,moms2) = rec_moments(v[(cut+1):n],ord)
		# now stitch them together
		
		n = n1 + n2  # redundant
		del21 = mu2 - mu1
		mu = mu1 + (n2 * del21 / n)   # eqn (II.3)
		sum_mom = moms1 + moms2
		nfoo = n1 * n2 * del21 / n
		for pp=2:ord
			sum_mom[pp-1] += (nfoo^pp) * ((n2^(1.0-pp)) - ((-n1)^(1.0-pp)))
			if (pp > 2)
				for kk=1:(pp - 2)
					sum_mom[pp-1] += binomial(pp,kk) * (del21^kk) * (((convert(Float64,-n2/n))^kk) * moms1[pp-kk-1] + 
																													 ((convert(Float64,n1/n))^kk) * moms2[pp-kk-1])
				end
			end
		end
		return (n, mu, sum_mom)
	end
end

function std3(v::AbstractVector)
	# return the standard deviation, the mean and the dof
	a,b,c = rec_moments(v,2)
	return (sqrt(c[1]/(a-1)),b,a)
end

function skew4(v::AbstractVector)
	# return the skew, the standard deviation, the mean, and the dof
	a,b,c = rec_moments(v,3)
	return (sqrt(a) * c[2] / (c[1]^1.5),
					sqrt(c[1]/(a-1)),b,a)
end

function kurt5(v::AbstractVector)
	# return the //excess// kurtosis, skew, standard deviation, mean, and the dof
	a,b,c = rec_moments(v,4)
	return ((a * c[3] / (c[1]^2.0)) - 3.0,
					sqrt(a) * c[2] / (c[1]^1.5),
					sqrt(c[1]/(a-1)),b,a)
end

function wel_std3(v::AbstractVector)
	# return the standard deviation, the mean and the dof
	# via Welford's method
	# see http://www.johndcook.com/standard_deviation.html
	df = numel(v)
  mu = v[1]
	sg = 0.0
	for iii=2:df
		dels = v[iii] - mu
		mu = mu + dels / iii
		# more in line with Bennett
		#sg = sg + ((iii - 1.0) / iii) * dels ^ 2.0
		sg = sg + dels * (v[iii] - mu)
	end
	return (sqrt(sg/(df-1)),mu,df)
end


# vim:ts=2:sw=2:tw=79:fdm=indent:cms=#%s:syntax=julia:filetype=julia:ai:si:cin:nu:fo=croql:cino=p0t0c5(0:

code for robust computation of moments in Julia via Welford/Bennett et. al/Terriberry method. My Julia foo is lacking, so this is rather slow.