The weighted mean¶

The weighted arithmetic mean of a table's members of the item's own, sum(w v) / sum(w), each sum in the fixed tree; an item with no member, or whose weights sum to nothing, reads 0.

Over a game's whole table of members (as lusine.finance.fair_value takes its loadings): each member's owner (the item it belongs to), its value and its weight. Equal weights give the plain mean; a game's linear or exponential weighting is the weights it writes. Generic: a value in any unit would be the mean in that unit, but a kernel's list fields are typed, so the value and the mean are dimensionless here and a game scales. Source: Bevington & Robinson, Data Reduction and Error Analysis for the Physical Sciences, ch. 4 (the weighted mean). Checked against a reference in numpy, which has its own tests.

Source

Bevington & Robinson, Data Reduction and Error Analysis for the Physical Sciences, ch. 4 (the weighted mean).

Proof

reference, invariant, limit.

weighted_mean.trait

namespace lusine.statistics;

float weighted_mean(list<struct { relation owner; float value; float weight; }> members)
{
  let total = sum(weight) where owner == it;
  if (total <= 0.0) return 0.0;
  return (sum(value * weight) where owner == it) / total;
}

The weighted arithmetic mean of a table's members of the item's own, sum(w v) / sum(w), each sum in Trait's fixed tree; 0 for an item with no member or with weights summing to nothing. The kernel takes a game's whole table of members, list<struct { relation owner; float value; float weight; }>, as lusine.finance.fair_value takes its loadings.

In [1]:

import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
    root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
import matplotlib.pyplot as plt
from lusine import trait_run
F32 = np.float32
U = float(np.finfo(F32).eps) / 2        # binary32's unit roundoff, 2^-24

import os, sys root = os.getcwd() while not os.path.isfile(os.path.join(root, 'docsite.json')): root = os.path.dirname(root) sys.path.insert(0, os.path.join(root, 'tools')) import numpy as np import matplotlib.pyplot as plt from lusine import trait_run F32 = np.float32 U = float(np.finfo(F32).eps) / 2 # binary32's unit roundoff, 2^-24

The kernel¶

The law as the library states it.

In [2]:

print(trait_run.source('statistics/weighted_mean.trait'))

print(trait_run.source('statistics/weighted_mean.trait'))

// lusine.statistics.weighted_mean: the weighted arithmetic mean of a table's members of the item's own, sum(w v) /
// sum(w), each sum in the fixed tree; an item with no member, or whose weights sum to nothing, reads 0. Over a game's
// whole table of members (as lusine.finance.fair_value takes its loadings): each member's owner (the item it belongs
// to), its value and its weight. Equal weights give the plain mean; a game's linear or exponential weighting is the
// weights it writes. Generic: a value in any unit would be the mean in that unit, but a kernel's list fields are
// typed, so the value and the mean are dimensionless here and a game scales.
// Source: Bevington & Robinson, Data Reduction and Error Analysis for the Physical Sciences, ch. 4 (the weighted mean).
// Checked against a reference in numpy, which has its own tests.

namespace lusine.statistics;

float weighted_mean(list<struct { relation owner; float value; float weight; }> members)
{
  let total = sum(weight) where owner == it;
  if (total <= 0.0) return 0.0;
  return (sum(value * weight) where owner == it) / total;
}

The program¶

A game's tables, written out: the items, each with its mean; and the members, each with its owner, its value and its weight, the values and weights the host's. means(values, weights) runs one step and gives each item's mean; MEMBERS says which item each member belongs to.

In [3]:

ITEMS = ['plain', 'linear', 'lone', 'empty', 'weightless', 'many']
MEMBERS = ['plain'] * 5 + ['linear'] * 6 + ['lone'] + ['weightless'] * 3 + ['many'] * 300


def program(kernel='lusine.statistics.weighted_mean', extra=''):
    items = ', '.join(f'{{ name = "{i}" }}' for i in ITEMS)
    members = ',\n'.join(f'  {{ owner = items["{i}"] }}' for i in MEMBERS)
    return f'''namespace stat;

lusine.version version = {{ {{ number = 1 }} }};

struct item {{ text name; float mean; }};
table<item> items(capacity = 8, key = name) = {{ {items} }};
struct member {{ relation<items> owner; host float value; host float weight; }};
table<member> members(capacity = {len(MEMBERS)}) = {{
{members},
}};
{extra}
void main() {{ items.mean = {kernel}(members); }}
'''

LAW = program()
print(LAW[:900], '...')


def means(values, weights, text=LAW):
    """{item: its mean} after one step, values and weights one a member in MEMBERS' order."""
    run = trait_run.Program(['statistics/weighted_mean.trait'], text, {})
    run.store('stat.members', 'value', np.asarray(values, dtype=float))
    run.store('stat.members', 'weight', np.asarray(weights, dtype=float))
    out = run.steps('stat.main', 1, [('stat.items', 'mean')])[0][0]
    return dict(zip(ITEMS, out[:len(ITEMS)]))


def tree(xs):
    if len(xs) == 0:
        return F32(0.0)
    if len(xs) == 1:
        return F32(xs[0])
    mid = (len(xs) + 1) // 2
    return F32(tree(xs[:mid]) + tree(xs[mid:]))


def model(values, weights):
    """Each item's mean as the kernel computes it: its members' w v and w each in the fixed tree, in slot order, divided."""
    out = {}
    for i in ITEMS:
        v = [F32(x) for x, o in zip(values, MEMBERS) if o == i]
        w = [F32(x) for x, o in zip(weights, MEMBERS) if o == i]
        total = tree(w)
        out[i] = F32(0.0) if total <= 0.0 else F32(tree([F32(a * b) for a, b in zip(v, w)]) / total)
    return out

ITEMS = ['plain', 'linear', 'lone', 'empty', 'weightless', 'many'] MEMBERS = ['plain'] * 5 + ['linear'] * 6 + ['lone'] + ['weightless'] * 3 + ['many'] * 300 def program(kernel='lusine.statistics.weighted_mean', extra=''): items = ', '.join(f'{{ name = "{i}" }}' for i in ITEMS) members = ',\n'.join(f' {{ owner = items["{i}"] }}' for i in MEMBERS) return f'''namespace stat; lusine.version version = {{ {{ number = 1 }} }}; struct item {{ text name; float mean; }}; table<item> items(capacity = 8, key = name) = {{ {items} }}; struct member {{ relation<items> owner; host float value; host float weight; }}; table<member> members(capacity = {len(MEMBERS)}) = {{ {members}, }}; {extra} void main() {{ items.mean = {kernel}(members); }} ''' LAW = program() print(LAW[:900], '...') def means(values, weights, text=LAW): """{item: its mean} after one step, values and weights one a member in MEMBERS' order.""" run = trait_run.Program(['statistics/weighted_mean.trait'], text, {}) run.store('stat.members', 'value', np.asarray(values, dtype=float)) run.store('stat.members', 'weight', np.asarray(weights, dtype=float)) out = run.steps('stat.main', 1, [('stat.items', 'mean')])[0][0] return dict(zip(ITEMS, out[:len(ITEMS)])) def tree(xs): if len(xs) == 0: return F32(0.0) if len(xs) == 1: return F32(xs[0]) mid = (len(xs) + 1) // 2 return F32(tree(xs[:mid]) + tree(xs[mid:])) def model(values, weights): """Each item's mean as the kernel computes it: its members' w v and w each in the fixed tree, in slot order, divided.""" out = {} for i in ITEMS: v = [F32(x) for x, o in zip(values, MEMBERS) if o == i] w = [F32(x) for x, o in zip(weights, MEMBERS) if o == i] total = tree(w) out[i] = F32(0.0) if total <= 0.0 else F32(tree([F32(a * b) for a, b in zip(v, w)]) / total) return out

namespace stat;

lusine.version version = { { number = 1 } };

struct item { text name; float mean; };
table<item> items(capacity = 8, key = name) = { { name = "plain" }, { name = "linear" }, { name = "lone" }, { name = "empty" }, { name = "weightless" }, { name = "many" } };
struct member { relation<items> owner; host float value; host float weight; };
table<member> members(capacity = 315) = {
  { owner = items["plain"] },
  { owner = items["plain"] },
  { owner = items["plain"] },
  { owner = items["plain"] },
  { owner = items["plain"] },
  { owner = items["linear"] },
  { owner = items["linear"] },
  { owner = items["linear"] },
  { owner = items["linear"] },
  { owner = items["linear"] },
  { owner = items["linear"] },
  { owner = items["lone"] },
  { owner = items["weightless"] },
  { owner = items["weightless"] },
  { owner = items["weightless"] },
  { owner = items["many"] },
  { ...

Reference: Bevington's weighted mean, and the float64 reference¶

Three measurements, 10.2 ± 0.1, 9.8 ± 0.2 and 10.5 ± 0.5, weighted by 1/sigma^2 (Bevington & Robinson, ch. 4), average to 10.132 as plain's first three members (the other two weightless); linear's six members with weights 1 to 6, and many's 300 random ones, are the float64 reference's within the folds' rounding.

In [4]:

from lusine.kernels.statistics.weighted_mean import weighted_mean as reference
rng = np.random.default_rng(420)
values = np.concatenate([[10.2, 9.8, 10.5, 0.0, 0.0], np.arange(1.0, 7.0), [7.0], [1.0, 2.0, 3.0], rng.normal(0.0, 1.0, 300)])
weights = np.concatenate([[100.0, 25.0, 4.0, 0.0, 0.0], np.arange(1.0, 7.0), [0.5], [0.0, 0.0, 0.0], rng.uniform(0.0, 2.0, 300)])
got = means(values, weights)
assert abs(got['plain'] - 10.132) < 0.001, got['plain']
for i in ('plain', 'linear', 'many'):
    v = [x for x, o in zip(values, MEMBERS) if o == i]
    w = [x for x, o in zip(weights, MEMBERS) if o == i]
    r = reference(F32(v).astype(float), F32(w).astype(float))
    assert abs(got[i] - r) <= 64 * U * abs(r) + 64 * U, (i, got[i], r)
print({i: round(float(got[i]), 6) for i in ITEMS})

from lusine.kernels.statistics.weighted_mean import weighted_mean as reference rng = np.random.default_rng(420) values = np.concatenate([[10.2, 9.8, 10.5, 0.0, 0.0], np.arange(1.0, 7.0), [7.0], [1.0, 2.0, 3.0], rng.normal(0.0, 1.0, 300)]) weights = np.concatenate([[100.0, 25.0, 4.0, 0.0, 0.0], np.arange(1.0, 7.0), [0.5], [0.0, 0.0, 0.0], rng.uniform(0.0, 2.0, 300)]) got = means(values, weights) assert abs(got['plain'] - 10.132) < 0.001, got['plain'] for i in ('plain', 'linear', 'many'): v = [x for x, o in zip(values, MEMBERS) if o == i] w = [x for x, o in zip(weights, MEMBERS) if o == i] r = reference(F32(v).astype(float), F32(w).astype(float)) assert abs(got[i] - r) <= 64 * U * abs(r) + 64 * U, (i, got[i], r) print({i: round(float(got[i]), 6) for i in ITEMS})

{'plain': 10.131783, 'linear': 4.333333, 'lone': 7.0, 'empty': 0.0, 'weightless': 0.0, 'many': 0.012963}

Reference: bit for bit, the folds in binary32¶

Every item's mean is the binary32 model's bit for bit: its members' w v and w each added in Trait's fixed tree in slot order, then divided.

In [5]:

want = model(values, weights)
assert all(F32(got[i]) == want[i] for i in ITEMS), {i: (got[i], want[i]) for i in ITEMS if F32(got[i]) != want[i]}
print('6 items bit for bit with the binary32 model')

want = model(values, weights) assert all(F32(got[i]) == want[i] for i in ITEMS), {i: (got[i], want[i]) for i in ITEMS if F32(got[i]) != want[i]} print('6 items bit for bit with the binary32 model')

6 items bit for bit with the binary32 model

Invariant: equal weights, and the weights' scale¶

With every weight 1, many is the plain mean of its 300 values within the folds' rounding; the weights tripled, every mean is unchanged within a rounding; and each mean lies among its values.

In [6]:

plain = means(values, np.ones_like(weights))
v300 = F32(values[-300:]).astype(float)
assert abs(plain['many'] - np.mean(v300)) <= 64 * U
tripled = means(values, 3.0 * weights)
assert all(abs(tripled[i] - got[i]) <= 8 * U * abs(got[i]) + 8 * U for i in ITEMS)
for i in ('plain', 'linear', 'many'):
    v = [x for x, o in zip(values, MEMBERS) if o == i]
    assert min(v) - 1e-6 <= got[i] <= max(v) + 1e-6
print('plain mean of many:', float(plain['many']), '; tripled weights leave', {i: round(float(tripled[i]), 6) for i in ('plain', 'linear')})

plain = means(values, np.ones_like(weights)) v300 = F32(values[-300:]).astype(float) assert abs(plain['many'] - np.mean(v300)) <= 64 * U tripled = means(values, 3.0 * weights) assert all(abs(tripled[i] - got[i]) <= 8 * U * abs(got[i]) + 8 * U for i in ITEMS) for i in ('plain', 'linear', 'many'): v = [x for x, o in zip(values, MEMBERS) if o == i] assert min(v) - 1e-6 <= got[i] <= max(v) + 1e-6 print('plain mean of many:', float(plain['many']), '; tripled weights leave', {i: round(float(tripled[i]), 6) for i in ('plain', 'linear')})

plain mean of many: 0.061385080218315125 ; tripled weights leave {'plain': 10.131783, 'linear': 4.333333}

Limit: no member, weights of nothing, one member¶

empty has no member and reads 0; weightless's three members weigh nothing and it reads 0; lone's one member is its own value.

In [7]:

assert got['empty'] == 0.0 and got['weightless'] == 0.0 and got['lone'] == 7.0, got
print(got['empty'], got['weightless'], got['lone'])

assert got['empty'] == 0.0 and got['weightless'] == 0.0 and got['lone'] == 7.0, got print(got['empty'], got['weightless'], got['lone'])

0.0 0.0 7.0

Planted fault¶

A wrong kernel over the same list, written here in Trait and bound in the library's place: the weighted sum undivided. The weights' scale catches it.

In [8]:

UNDIVIDED = program('undivided', '''
float undivided(list<struct { relation owner; float value; float weight; }> members)
{
  return sum(value * weight) where owner == it;
}
''')
W = means(values, 3.0 * weights, text=UNDIVIDED)
assert not all(abs(W[i] - got[i]) <= 8 * U * abs(got[i]) + 8 * U for i in ITEMS)
print('PLANTED: the undivided sum is caught')

UNDIVIDED = program('undivided', ''' float undivided(list<struct { relation owner; float value; float weight; }> members) { return sum(value * weight) where owner == it; } ''') W = means(values, 3.0 * weights, text=UNDIVIDED) assert not all(abs(W[i] - got[i]) <= 8 * U * abs(got[i]) + 8 * U for i in ITEMS) print('PLANTED: the undivided sum is caught')

PLANTED: the undivided sum is caught