Why cmomy?

\[\newcommand{\ave}[1]{\langle{#1}\rangle}\]

We wish to calculate moments of the form \(\ave{x}\) and \(\ave{(\delta x)^k}\), \(\delta x = (x - \ave{x})\), where, in general we define averages of the form.

\[ \ave{x} = \frac{\sum_i w_i x_i}{\sum_i w_i} \]

\[ \ave{(\delta x)^k} = \frac{\sum_i w_i (x_i - \ave{x})^k}{\sum_i w_i} \]

where \(w_i\) are weights associated with each sample \(x_i\). In the simple case, \(w_i = 1\). To calculate a central moment with \(k > 1\), we have some choices.

1. First calculate \(\ave{x}\), then use the above formula.

• Downside: requires two passes through the data.

• Upside: numerically stable.

2. Expand out the central moment to raw moments. For example \(\ave{(x - \ave{x})^2} = \ave{x^n} - \ave{x}^2\). These raw moments can be calculated directly in one pass.

• Downside: numerically unstable (see here for some info).

• Upside: single pass for each of \(\ave{x^k}\).

An example

In what follows, we’ll arrange data in the from output[moment] where output[0] = {total weight}, output[1] = {average value}, output[k > 1] = {kth order central moment}. We can generate the central moments using the following functions

def cmom_1(x, axis=0, mom=3):
    """Create central moments using stable method"""
    # output shape
    shape = x.shape[:axis] + x.shape[axis + 1 :] + (mom + 1,)
    output = np.zeros(shape, dtype=x.dtype)

    output[..., 0] = x.shape[axis]
    output[..., 1] = x.mean(axis=axis)

    # moments [2 -> mom]
    output[..., 2:] = (
        (x - x.mean(axis=axis, keepdims=True))[..., None] ** np.arange(2, mom + 1)
    ).mean(axis=axis)
    return output


def cmom_2(x, axis=0, mom=3):
    if mom > 3:
        raise ValueError

    shape = x.shape[:axis] + x.shape[axis + 1 :] + (mom + 1,)

    output = np.zeros(shape, dtype=x.dtype)

    output[..., 0] = x.shape[axis]

    # higher order raw moments
    raws = (x[..., None] ** np.arange(0, mom + 1)).mean(axis=axis)

    output[..., 1] = raws[..., 1]
    if mom > 1:
        # <(x - <x>)**2> = <x**2> - <x>**2
        output[..., 2] = raws[..., 2] - raws[..., 1] ** 2

    if mom > 2:
        # <(x - <x>)**3> = <x**3> - 3<x**2><x> + 2<x>**3
        output[..., 3] = (
            raws[..., 3] - 3 * raws[..., 2] * raws[..., 1] + 2 * raws[..., 1] ** 3
        )

    return output

import numpy as np

rng = np.random.default_rng(seed=0)

n = 1000
x = rng.random(n)

m1 = cmom_1(x)
m2 = cmom_2(x)
# same values
np.testing.assert_allclose(m1, m2)

OK, great. But what if we have unscaled data? For example, setting \(y = a x + b\) gives \(\ave{y} = a \ave{x} + b\), but for central moments, \(\ave{(y - \ave{y})^n} = \ave{(a x - b - \ave{a x - b})^n} = a^n \ave{(\delta x)^n}\)

# large value, small fluctuation
a = 0.1
b = 1e4
y = a * x + b

# expected value after shift
expected = m1.copy()
expected[..., 1] = a * expected[..., 1] + b
expected[..., 2:] *= a ** np.arange(2, 4)


m1_shift = cmom_1(y)
m2_shift = cmom_2(y)

m1_error = m1_shift - expected
m2_error = m2_shift - expected

print(
    f"""
abs error using central moments: {m1_error}
rel error using central moments: {m1_error / expected}
abs error using raw moments    : {m2_error}
rel error using raw moments    : {m2_error / expected}
"""
)

abs error using central moments: [ 0.0000e+00 -1.8190e-12 -8.8373e-16  3.5155e-15]
rel error using central moments: [ 0.0000e+00 -1.8190e-16 -1.0912e-12 -1.6856e-09]
abs error using raw moments    : [ 0.0000e+00  7.2760e-12 -7.2259e-08  9.7865e-04]
rel error using raw moments    : [ 0.0000e+00  7.2759e-16 -8.9221e-05 -4.6925e+02]

So what?

OK, so central moments have some advantages. But the two pass thing makes life difficult. For example, what if we are running a slow experiment or simulation? We could collect all the samples for x, then do the two pass algorithm at the end. But it would be nice to calculate averages on the fly. This is why we often turn to raw moments. There easy to accumulate as data comes in, while central moments, at least as formulated above, are not.

Thankfully, smart people have figured out how to calculate central moments in a one-pass way. Moreover, these methods allow central moments to be combined. This makes resampling from central moments possible. We don’t go into detail on these formulas, but point the interested reader to the this article, and citations within.

The cmomy package has a variety of routines to work with central moments. This include

• Create central moments accumulator object in a variety of ways.

• add data to object

• resample data

• extract central and raw moments

The primary object for doing all this is CentralMoments.

import cmomy

cx = cmomy.CentralMoments.from_vals(x, axis=0, mom=3)
cy = cmomy.CentralMoments.from_vals(y, axis=0, mom=3)


cx

array([ 1.0000e+03,  5.1691e-01,  8.0989e-02, -2.0856e-03])

CentralMoments

• mom: 3

• Info: (0)

• 1e+03 0.5169 0.08099 -0.002086

  array([ 1.000000e+03,  5.169063e-01,  8.098908e-02, -2.085567e-03])

You can access the underlying data using the cmomy.CentralMoments.values or cmomy.CentralMoments.data attributes.

print(cx.values)

[ 1.0000e+03  5.1691e-01  8.0989e-02 -2.0856e-03]

Basically, the CentralMoments object is just a collection of routines around data, which is an array of central moments (with the convention that data[..., 0] is the total weight, data[..., 1] is the average, and data[..., k > 1] is the kth central moment). It gives the same results as calculating central moments directly:

np.testing.assert_allclose(cx.values, m1)
np.testing.assert_allclose(cy.values, m1_shift)

Nice! But what can that’s not all that special. The real power comes in for multidimensional data. Suppose that we have separate samples of the data \(\{x_a, x_b, ... \}\), and want to get there central moments. We can do the following.

x = rng.random((100, 5))
m1 = cmom_1(x, axis=0, mom=4)
m1

array([[ 1.0000e+02,  4.7703e-01,  8.6779e-02, -2.0776e-03,  1.3134e-02],
       [ 1.0000e+02,  4.7149e-01,  8.9923e-02,  4.1844e-06,  1.3357e-02],
       [ 1.0000e+02,  4.9054e-01,  8.3374e-02,  6.0526e-03,  1.2658e-02],
       [ 1.0000e+02,  5.0095e-01,  9.1982e-02, -1.1610e-03,  1.4399e-02],
       [ 1.0000e+02,  5.1960e-01,  9.1251e-02, -2.4683e-03,  1.4763e-02]])

Wrapping the data in a central moments object gives

c = cmomy.CentralMoments.from_data(m1, mom_ndim=1)

# reduce the data along a dimension to get

print(c.reduce(axis=0).values)

# verify that this is the same as just accumulating all the data
print(cmom_1(x.reshape(-1), mom=4))

[5.0000e+02 4.9192e-01 8.8959e-02 1.3873e-04 1.3778e-02]
[5.0000e+02 4.9192e-01 8.8959e-02 1.3873e-04 1.3778e-02]

Excellent! We can use CentralMoments to accumulate data on the go

c = cmomy.CentralMoments.zeros(mom=4)
for i in range(5):
    # note: pushing multiple values
    c.push_vals(x[:, i])
c

array([5.0000e+02, 4.9192e-01, 8.8959e-02, 1.3873e-04, 1.3778e-02])

CentralMoments

• mom: 4

• Info: (0)

• 500.0 0.4919 0.08896 0.0001387 0.01378

  array([5.000000e+02, 4.919195e-01, 8.895944e-02, 1.387264e-04, 1.377840e-02])

Or, you can push all the values one at a time!

c.zero()
for xx in x.reshape(-1):
    # note: pushing single value
    c.push_val(xx)

c

array([5.0000e+02, 4.9192e-01, 8.8959e-02, 1.3873e-04, 1.3778e-02])

CentralMoments

• mom: 4

• Info: (0)

• 500.0 0.4919 0.08896 0.0001387 0.01378

  array([5.000000e+02, 4.919195e-01, 8.895944e-02, 1.387264e-04, 1.377840e-02])

The values can even be differently weighted.

xs = np.split(x.reshape(-1), [50, 150, 300])

c.zero()
for v in xs:
    c.push_vals(v)
c

array([5.0000e+02, 4.9192e-01, 8.8959e-02, 1.3873e-04, 1.3778e-02])

CentralMoments

• mom: 4

• Info: (0)

• 500.0 0.4919 0.08896 0.0001387 0.01378

  array([5.000000e+02, 4.919195e-01, 8.895944e-02, 1.387264e-04, 1.377840e-02])

Alternatively, you can add objects together.

cs = [cmomy.CentralMoments.from_vals(v, mom=4) for v in xs]
print(cs)

c_tot = cs[0] + cs[1] + cs[2] + cs[3]
c_tot

[<CentralMoments(val_shape=(), mom=(4,))>
array([5.0000e+01, 4.1743e-01, 7.3610e-02, 4.5490e-03, 1.0645e-02]), <CentralMoments(val_shape=(), mom=(4,))>
array([ 1.0000e+02,  5.1217e-01,  9.0418e-02, -3.0034e-03,  1.3257e-02]), <CentralMoments(val_shape=(), mom=(4,))>
array([ 1.5000e+02,  4.8912e-01,  8.3731e-02, -6.4426e-04,  1.2724e-02]), <CentralMoments(val_shape=(), mom=(4,))>
array([2.0000e+02, 5.0252e-01, 9.4279e-02, 1.8869e-04, 1.5262e-02])]

array([5.0000e+02, 4.9192e-01, 8.8959e-02, 1.3873e-04, 1.3778e-02])

CentralMoments

• mom: 4

• Info: (0)

• 500.0 0.4919 0.08896 0.0001387 0.01378

  array([5.000000e+02, 4.919195e-01, 8.895944e-02, 1.387264e-04, 1.377840e-02])

Or more simply:

print(sum(cs, start=cs[0].zeros_like()).values)

[5.0000e+02 4.9192e-01 8.8959e-02 1.3873e-04 1.3778e-02]

Very cool!

weighted averages

cmomy is setup to work with arbitrary weights. We saw this work above when we considered unequal sample sizes. For example,

def get_cmom_with_weights(x, w, axis=0, mom=3):
    """Create central moments using stable method"""
    # output shape
    shape = x.shape[:axis] + x.shape[axis + 1 :] + (mom + 1,)
    output = np.zeros(shape, dtype=x.dtype)

    w_sum = w.sum(axis=axis)
    w_norm = 1.0 / w_sum

    output[..., 0] = w_sum
    output[..., 1] = ((w * x).sum(axis=axis)) * w_norm

    # moments [2 -> mom]
    output[..., 2:] = (
        w[..., None]
        * (x - (w * x).sum(axis=axis, keepdims=True) * w_norm)[..., None]
        ** np.arange(2, mom + 1)
    ).sum(axis=axis) * w_norm

    return output

x = rng.random(100)
w = rng.random(100)

moments = get_cmom_with_weights(x, w)
moments

cmomy.CentralMoments.from_vals(x, w=w, mom=3)

array([4.8363e+01, 5.0295e-01, 9.2209e-02, 6.3664e-04])

CentralMoments

• mom: 3

• Info: (0)

• 48.36 0.503 0.09221 0.0006366

  array([4.836319e+01, 5.029529e-01, 9.220858e-02, 6.366441e-04])

comoments

cmomy can also handle comoments (covariance, etc), like the \(\ave{(\delta x)^i (\delta y)^j}\)
For example

y = rng.random(100)

cmomy.CentralMoments.from_vals(x=(x, y), w=w, mom=(2, 2))

array([[ 4.8363e+01,  4.8516e-01,  7.5856e-02],
       [ 5.0295e-01, -1.2237e-02,  3.2919e-04],
       [ 9.2209e-02, -4.2304e-03,  6.5653e-03]])

CentralMoments

• mom: (2, 2)

• Info: (0)

• 48.36 0.4852 0.07586 0.503 ... 0.0003292 0.09221 -0.00423 0.006565

  array([[ 4.836319e+01,  4.851622e-01,  7.585550e-02],
         [ 5.029529e-01, -1.223740e-02,  3.291919e-04],
         [ 9.220858e-02, -4.230387e-03,  6.565301e-03]])

All the special sauce works for central moments or comoments.

xarray DataArray support

cmomy also works with xarray.DataArray objects. For example

import xarray as xr

x = xr.DataArray(rng.random((100, 2, 3)), dims=["rec", "a", "b"])
x.head(3)  # type: ignore[attr-defined]

<xarray.DataArray (rec: 3, a: 2, b: 3)> Size: 144B
array([[[0.0625, 0.5585, 0.0184],
        [0.4284, 0.4153, 0.0854]],

       [[0.0821, 0.587 , 0.0041],
        [0.7317, 0.3686, 0.2628]],

       [[0.951 , 0.0224, 0.627 ],
        [0.0179, 0.382 , 0.3126]]])
Dimensions without coordinates: rec, a, b

xarray.DataArray

• rec: 3
• a: 2
• b: 3

• 0.06247 0.5585 0.01837 0.4284 0.4153 ... 0.627 0.01787 0.382 0.3126

  array([[[0.0625, 0.5585, 0.0184],
          [0.4284, 0.4153, 0.0854]],
  
         [[0.0821, 0.587 , 0.0041],
          [0.7317, 0.3686, 0.2628]],
  
         [[0.951 , 0.0224, 0.627 ],
          [0.0179, 0.382 , 0.3126]]])

• Coordinates: (0)
• Indexes: (0)
• Attributes: (0)

To create a cmomy object wrapping a xarray.DataArray, use xCentralMoments

c = cmomy.xCentralMoments.from_vals(x, dim="rec", mom=3)
c

<xarray.DataArray (a: 2, b: 3, mom_0: 4)> Size: 192B
array([[[1.0000e+02, 4.0657e-01, 9.0470e-02, 1.3814e-02],
        [1.0000e+02, 4.9516e-01, 7.3811e-02, 1.1677e-03],
        [1.0000e+02, 4.6592e-01, 8.1210e-02, 3.3827e-03]],

       [[1.0000e+02, 4.8831e-01, 8.2432e-02, 1.6465e-04],
        [1.0000e+02, 4.6920e-01, 9.2950e-02, 5.5999e-03],
        [1.0000e+02, 4.8891e-01, 7.6961e-02, 2.3232e-03]]])
Dimensions without coordinates: a, b, mom_0

xCentralMoments

• mom: 3
• val_shape: (2, 3)

• Info: (2)

  mom_dims :

  mom_0

  val_dims :

  (a, b)

<xarray.DataArray (a: 2, b: 3, mom_0: 4)> Size: 192B
array([[[1.0000e+02, 4.0657e-01, 9.0470e-02, 1.3814e-02],
        [1.0000e+02, 4.9516e-01, 7.3811e-02, 1.1677e-03],
        [1.0000e+02, 4.6592e-01, 8.1210e-02, 3.3827e-03]],

       [[1.0000e+02, 4.8831e-01, 8.2432e-02, 1.6465e-04],
        [1.0000e+02, 4.6920e-01, 9.2950e-02, 5.5999e-03],
        [1.0000e+02, 4.8891e-01, 7.6961e-02, 2.3232e-03]]])
Dimensions without coordinates: a, b, mom_0

xarray.DataArray

• a: 2
• b: 3
• mom_0: 4

• 100.0 0.4066 0.09047 0.01381 100.0 ... 100.0 0.4889 0.07696 0.002323

  array([[[1.0000e+02, 4.0657e-01, 9.0470e-02, 1.3814e-02],
          [1.0000e+02, 4.9516e-01, 7.3811e-02, 1.1677e-03],
          [1.0000e+02, 4.6592e-01, 8.1210e-02, 3.3827e-03]],
  
         [[1.0000e+02, 4.8831e-01, 8.2432e-02, 1.6465e-04],
          [1.0000e+02, 4.6920e-01, 9.2950e-02, 5.5999e-03],
          [1.0000e+02, 4.8891e-01, 7.6961e-02, 2.3232e-03]]])

• Coordinates: (0)
• Indexes: (0)
• Attributes: (0)

Everything that CentralMoments can do, xCentralMoments can do.