|
|
|
@ -1,6 +1,9 @@ |
|
|
|
|
# $Id: avg.py,v 1.1 2010-10-07 15:56:55 wirawan Exp $ |
|
|
|
|
# Create: 20090112 |
|
|
|
|
|
|
|
|
|
# TODO: Use less roundoff-error algorithm found in: |
|
|
|
|
# http://en.wikipedia.org/wiki/Standard_deviation#Weighted_calculation |
|
|
|
|
|
|
|
|
|
import math |
|
|
|
|
|
|
|
|
|
class avg(object): |
|
|
|
@ -27,3 +30,48 @@ class avg(object): |
|
|
|
|
# Operator overloading: |
|
|
|
|
avg.__iadd__ = avg.add |
|
|
|
|
avg.__call__ = avg.mean |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
class weighted_stats(object): |
|
|
|
|
def __init__(self, a=None, weights=None): |
|
|
|
|
if a != None: |
|
|
|
|
get_sample_stats(a, weights, result=self) |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
def get_sample_stats(a, weights=None, result=None): |
|
|
|
|
"""Applies basic statistics (average, variance, standard deviation) |
|
|
|
|
for a given sample, optionally with a weight. |
|
|
|
|
""" |
|
|
|
|
from numpy import asarray, count_nonzero, product, sqrt, sum, nan |
|
|
|
|
a = asarray(a) |
|
|
|
|
if result == None: |
|
|
|
|
r = weighted_stats() |
|
|
|
|
else: |
|
|
|
|
r = result |
|
|
|
|
if weights == None: |
|
|
|
|
r.s0 = product(a.shape) |
|
|
|
|
r.s1 = sum(a) |
|
|
|
|
r.s2 = sum(a**2) |
|
|
|
|
r.N = r.s0 |
|
|
|
|
r.N_nz = r.s0 |
|
|
|
|
else: |
|
|
|
|
weights = asarray(weights) |
|
|
|
|
r.s0 = sum(weights) |
|
|
|
|
r.s1 = sum(weights * a) |
|
|
|
|
r.s2 = sum(weights * a**2) |
|
|
|
|
r.N_nz = count_nonzero(weights) |
|
|
|
|
r.N = product(a.shape) |
|
|
|
|
# Ref: http://en.wikipedia.org/wiki/Standard_deviation#Weighted_calculation |
|
|
|
|
# TODO: Use less roundoff-error algorithm found therein. |
|
|
|
|
r.avg = r.s1 / r.s0 |
|
|
|
|
r.var_pop = (r.s0 * r.s2 - r.s1**2) / (r.s0**2) |
|
|
|
|
r.std_pop = sqrt(r.var_pop) |
|
|
|
|
if r.N_nz > 1: |
|
|
|
|
r.var_samp = (r.N_nz / (r.N_nz - 1.0)) * r.var_pop |
|
|
|
|
r.std_samp = sqrt(r.var_samp) |
|
|
|
|
else: |
|
|
|
|
r.var_samp = nan |
|
|
|
|
r.std_samp = nan |
|
|
|
|
|
|
|
|
|
return r |
|
|
|
|
|
|
|
|
|