Sensor placement for wearable fall detection, under honest validation¶

CMS7503 Machine Learning — Assignment 2 (Development Task)

This notebook is the evidence base for the accompanying report. It re-runs the sensor-placement comparison on the UCI Simulated Falls and Daily Living Activities dataset that Assignment 1 committed to, with one change that decides everything: every accuracy here is measured under leave-one-subject-out (LOSO) cross-validation, so each figure describes performance on a person the model has never seen. The published results on this dataset (Özdemir & Barshan, 2014; Özdemir, 2016; Ntanasis et al., 2017) all used k-fold cross-validation, which lets repetitions from the same volunteer sit in both the training and test folds. A deployed device always meets a new wearer, so that is the test I run.

The question the notebook exists to answer is narrow and was set in Assignment 1: under subject-independent validation, does the wrist come within about two to three points of the best trunk or thigh placement? If it does, the everyday acceptability of a wristband may outweigh a small accuracy cost, and the best placement for a real product diverges from the best placement on paper. If the gap stays wide, the case for a trunk or thigh device holds.

The notebook runs top to bottom and regenerates every number and figure the report cites.

In [1]:
import os, re, glob, time, warnings
import numpy as np
import pandas as pd
from scipy import stats, fft
from scipy.stats import wilcoxon, t as tdist
import matplotlib.pyplot as plt
import matplotlib as mpl

warnings.filterwarnings("ignore")
np.random.seed(0)

DATA = "data/Tests"          # extracted 17-subject release (see Section 1)
CACHE = "data/cache.npz"     # feature matrices, rebuilt below if absent

# a small, colour-blind-safe styling so the figures read as one set
mpl.rcParams.update({
    "figure.dpi": 110, "axes.spines.top": False, "axes.spines.right": False,
    "axes.grid": True, "grid.alpha": 0.25, "font.size": 11,
})
LOC_ORDER = ["head", "chest", "waist", "wrist", "thigh", "ankle"]
PALETTE = {"head": "#7b6cf6", "chest": "#4c9be8", "waist": "#2bb2a8",
           "wrist": "#e8734c", "thigh": "#6bbf59", "ankle": "#b58a3c"}

1. The data, and an audit before any modelling¶

The dataset holds 17 volunteers, each performing 20 fall types and 16 activities of daily living (ADLs) with roughly five repetitions, while wearing six Xsens MTw units on the head, chest, waist, right wrist, right thigh and right ankle. I load the raw Tests archive rather than the convenient ucimlrepo one-liner on purpose: that one-liner returns a single flat 138-feature table with no subject identifiers and no per-sensor split, and this project cannot be done without either. LOSO needs the subject label; a placement comparison needs the sensor separated out. So I parse the raw signal files, which carry both.

Two things must be confirmed before anything else. First, that this really is the 17-subject release — Assignment 1's framing depends on it, and the earlier papers worked from a smaller 14-subject version. Second, the class balance and the presence of missing values, both of which change how the rest of the notebook has to behave.

In [2]:
# The dataset documents a fixed sensor-serial -> body-location map (Information.pdf),
# and encodes the activity in the folder name: 8xx = ADL, 9xx = fall.
SENSOR_LOCATION = {"340506": "head", "340527": "chest", "340535": "waist",
                   "340537": "wrist", "340539": "thigh", "340540": "ankle"}
FS = 25.0                       # Hz, stated in every file header ("Update Rate: 25.0Hz")
USE_COLS = ["Acc_X", "Acc_Y", "Acc_Z", "Gyr_X", "Gyr_Y", "Gyr_Z"]

def read_sensor_file(path):
    # four '//' comment lines, then a tab-separated table with a trailing tab
    df = pd.read_csv(path, sep="\t", skiprows=4, engine="c")
    return df.loc[:, ~df.columns.str.startswith("Unnamed")]

def trial_index(root=DATA):
    rows = []
    for subj in sorted(os.listdir(root)):
        base = os.path.join(root, subj, "Testler Export")
        if not (subj.isdigit() and os.path.isdir(base)):
            continue
        for act_dir in sorted(os.listdir(base)):
            m = re.match(r"(\d{3})-", act_dir)
            if not m:
                continue
            act = int(m.group(1))
            for test_dir in sorted(glob.glob(os.path.join(base, act_dir, "Test_*"))):
                for f in glob.glob(os.path.join(test_dir, "*.txt")):
                    serial = os.path.splitext(os.path.basename(f))[0]
                    if serial in SENSOR_LOCATION:
                        rows.append(dict(subject=int(subj), activity=act,
                                         test=os.path.basename(test_dir),
                                         location=SENSOR_LOCATION[serial],
                                         path=f, is_fall=int(act >= 900)))
    return pd.DataFrame(rows)

index = trial_index()
n_subjects = index.subject.nunique()
print("subjects:", sorted(index.subject.unique()))
print("n_subjects =", n_subjects)
assert n_subjects == 17, "Not the 17-subject release — stop; A1's framing assumes 17."
print("sensor locations:", sorted(index.location.unique()))
print("rows (one per subject x activity x test x sensor):", len(index))
subjects: [np.int64(101), np.int64(102), np.int64(103), np.int64(104), np.int64(105), np.int64(106), np.int64(107), np.int64(108), np.int64(109), np.int64(110), np.int64(203), np.int64(204), np.int64(205), np.int64(206), np.int64(207), np.int64(208), np.int64(209)]
n_subjects = 17
sensor locations: ['ankle', 'chest', 'head', 'thigh', 'waist', 'wrist']
rows (one per subject x activity x test x sensor): 19950

The assertion is deliberate: if the archive were the 14-subject release the notebook stops here rather than quietly reporting on a different sample than the report claims.

In [3]:
# Class balance, at the level of one trial (all six sensors of one repetition share a label)
trials = index.drop_duplicates(["subject", "activity", "test"])
counts = trials.is_fall.value_counts().rename({0: "ADL", 1: "fall"})
print(counts)
print("fall fraction: %.3f" % trials.is_fall.mean())
is_fall
fall    1843
ADL     1483
Name: count, dtype: int64
fall fraction: 0.554

Falls are the majority class here, at roughly 55%, because the protocol has 20 fall types against 16 ADLs. That is the reverse of reality, where falls are rare, and it is the first reason accuracy alone is a poor headline: a classifier is rewarded for the common class, and here the common class is the fall. I report per-class metrics throughout for this reason, and I keep the natural balance rather than resampling, so the numbers stay comparable across the six locations. Recall on the fall class is sensitivity; I keep sklearn's column name recall in the tables and use the clinical term in the discussion, since Assignment 1 framed the problem in those terms.

The dataset documentation gives 3,060 trials; parsing the archive yields 3,326. The difference is real, not a parsing error: repetitions are not a flat five per activity. Most subject-activity pairs have five recordings and many have six, with a handful of seven or eight, for an average of about 5.4 — the documentation's round figure assumes five throughout and so undercounts. I report the number I actually parse and note the discrepancy rather than reconciling it silently.

In [4]:
# Missing values: the UCI documentation says "Rarely". Quantify it across every file
# so the handling choice is evidence-based, not assumed.
miss_files, miss_cells, checked = 0, 0, 0
empty = 0
sample = index.sample(frac=1.0, random_state=0)   # every file, shuffled
for p in sample.path:
    try:
        d = read_sensor_file(p)
    except Exception:
        empty += 1; continue
    if len(d) == 0:
        empty += 1; continue
    if set(USE_COLS).issubset(d.columns):
        checked += 1
        na = d[USE_COLS].isna().to_numpy().sum()
        if na:
            miss_files += 1; miss_cells += int(na)
print(f"files checked: {checked}")
print(f"files with >=1 missing signal value: {miss_files}")
print(f"total missing signal cells: {miss_cells}")
print(f"empty / unreadable files: {empty}  ({empty/len(index)*100:.2f}% of all files)")
files checked: 19928
files with >=1 missing signal value: 2765
total missing signal cells: 45606
empty / unreadable files: 22  (0.11% of all files)

Missing values are rare at the level of an individual sample — a fraction of a percent of all readings — though they are scattered across a meaningful share of files, and a small number of files are empty or truncated. Both are handled at the point of feature extraction (Section 3): missing samples are linearly interpolated within a channel and any residual filled at the ends, and the empty files are dropped and counted rather than allowed to poison a window. Every retained trial still carries its subject identifier and its sensor location — the two hinges — because those come from the folder path, not the signal.

2. From a trial to a labelled feature vector¶

This is the Activity Recognition Chain (Bulling et al., 2014) applied concretely: segment each recording, extract features, attach the label, the subject and the sensor. Two design choices matter and I make them explicitly.

Window length. Özdemir & Barshan used a four-second window centred on the peak total acceleration. At 25 Hz that is 100 samples. I mirror it: it is long enough to contain the impact and the settling that follows, and short enough that a single ADL or fall dominates the window.

Where to centre the window. I centre each sensor on its own peak acceleration magnitude, not on a peak shared across the body. A wrist and a waist do not reach their largest acceleration at the same instant during a fall — the limb whips while the trunk drops — so a shared centre would clip the very part of the distal signal that carries the information. Centring each sensor on its own peak is also what makes the six locations directly comparable, since each is given its best four seconds.

One limitation follows from this choice and belongs here rather than buried in the conclusion. Centring on the peak assumes the fall has already been located in the stream, which a continuously running device cannot do; a deployed system has to detect a candidate event first and classify it second. Every figure below therefore compares placements under identical and favourable segmentation. That is the right control for the question I am asking — which location carries the most information — and it is an upper bound on what a live system would achieve. Assignment 1 identified this ceiling in Özdemir and Barshan's protocol, and I inherit it deliberately to keep the comparison commensurable with the published results.

In [5]:
WINDOW_S = 4.0
WIN = int(WINDOW_S * FS)         # 100 samples

def window_on_peak(sig):
    # sig: (N, 6) = acc xyz, gyr xyz. Centre a WIN-sample window on peak acc magnitude.
    accmag = np.sqrt((sig[:, 0:3] ** 2).sum(1))
    peak = int(np.argmax(accmag))
    half = WIN // 2
    lo, hi = peak - half, peak - half + WIN
    if lo < 0:            lo, hi = 0, WIN
    if hi > len(sig):     lo, hi = len(sig) - WIN, len(sig)
    return sig[max(lo, 0):hi]

3. Features that can be explained¶

Assignment 1 argued that with only 17 subjects the appropriate choice is classical, feature-based machine learning rather than a data-hungry deep network. That argument only holds if the features are ones I can actually reason about, so I keep them standard and interpretable. From each of the six signal channels (three acceleration axes, three gyroscope axes) I take eight time-domain descriptors and three frequency-domain descriptors, and I add the same set for the acceleration signal-magnitude vector, which is orientation-independent and captures the overall intensity of movement.

Domain Features
Time mean, standard deviation, min, max, RMS, mean absolute deviation, skewness, kurtosis
Frequency dominant frequency, spectral energy, spectral entropy

That is 11 features per channel across seven channels (six sensor axes plus the magnitude vector), 77 in total. I exclude the magnetometer: it measures orientation relative to the earth's field, drifts indoors near metal and furniture, and carries little that distinguishes a fall from an ADL — including it would add noise I would then have to explain away.

In [6]:
def channel_features(x, fs=FS):
    f = {}
    f["mean"] = x.mean(); f["std"] = x.std(); f["min"] = x.min(); f["max"] = x.max()
    f["rms"] = np.sqrt((x ** 2).mean())
    f["mad"] = np.mean(np.abs(x - x.mean()))
    f["skew"] = stats.skew(x); f["kurt"] = stats.kurtosis(x)
    spec = np.abs(fft.rfft(x - x.mean()))
    freqs = fft.rfftfreq(len(x), d=1.0 / fs)
    power = spec ** 2
    f["dom_freq"] = freqs[np.argmax(spec)] if len(spec) else 0.0
    f["spec_energy"] = power.sum()
    pn = power / power.sum() if power.sum() > 0 else power
    f["spec_entropy"] = float(-(pn[pn > 0] * np.log2(pn[pn > 0])).sum())
    return f

def trial_features(win):
    feats = {}
    for i, name in enumerate(USE_COLS):
        for k, v in channel_features(win[:, i]).items():
            feats[f"{name}_{k}"] = v
    smv = np.sqrt((win[:, 0:3] ** 2).sum(1))
    for k, v in channel_features(smv).items():
        feats[f"smv_{k}"] = v
    return feats

Building the six feature matrices¶

The loop below turns every trial into a feature vector, once per sensor location, carrying the binary label and the subject identifier alongside. Building all six matrices reads roughly twenty thousand files and takes a few minutes, so the result is cached; the full build code is here and runs whenever the cache is absent, so the notebook stays reproducible from the raw data. Missing values are interpolated inside each channel and empty files are dropped and counted, exactly as described above.

In [7]:
def build_location_matrix(index, location):
    # returns both the feature matrix and the raw windows (the CNN in Section 8 needs
    # the latter); every trial is padded to >= WIN before windowing, so window_on_peak
    # always yields exactly WIN rows and np.array(W) is a clean (n, WIN, 6) array.
    sub = index[index.location == location]
    X, y, g, W, dropped = [], [], [], [], 0
    for _, r in sub.iterrows():
        df = read_sensor_file(r.path)
        if not set(USE_COLS).issubset(df.columns) or len(df) == 0:
            dropped += 1; continue
        sig = df[USE_COLS].to_numpy(dtype=float)
        if np.isnan(sig).any():
            sig = pd.DataFrame(sig).interpolate(limit_direction="both").to_numpy()
            sig = np.nan_to_num(sig, nan=0.0)
        if len(sig) < WIN:
            sig = np.pad(sig, ((0, WIN - len(sig)), (0, 0)), mode="edge")
        win = window_on_peak(sig)                 # (WIN, 6), same window feeds both paths
        X.append(trial_features(win)); W.append(win)
        y.append(r.is_fall); g.append(r.subject)
    return pd.DataFrame(X), np.array(y), np.array(g), np.array(W), dropped

if os.path.exists(CACHE):
    d = np.load(CACHE, allow_pickle=True)
    feat_names = list(d["feat_names"])
    feats = {loc: (d[f"X_{loc}"], d[f"y_{loc}"], d[f"g_{loc}"]) for loc in LOC_ORDER}
    raws = {loc: (d[f"W_{loc}"], d[f"yW_{loc}"], d[f"gW_{loc}"]) for loc in LOC_ORDER}
    print("loaded cached features")
else:
    feats, raws, feat_names, blob = {}, {}, None, {}
    for loc in LOC_ORDER:
        X, y, g, W, dropped = build_location_matrix(index, loc)
        feat_names = feat_names or list(X.columns)
        feats[loc] = (X.to_numpy(), y, g)
        raws[loc] = (W, y, g)
        blob[f"X_{loc}"], blob[f"y_{loc}"], blob[f"g_{loc}"] = X.to_numpy(), y, g
        blob[f"W_{loc}"], blob[f"yW_{loc}"], blob[f"gW_{loc}"] = W, y, g
        print(f"{loc}: X={X.shape}  raw={W.shape}  dropped={dropped}")
    os.makedirs(os.path.dirname(CACHE), exist_ok=True)
    np.savez_compressed(CACHE, feat_names=np.array(feat_names), **blob)
    print("cache written to", CACHE)

for loc in LOC_ORDER:
    X, y, g = feats[loc]
    print(f"{loc:6s}  X={X.shape}  fall_rate={y.mean():.3f}  subjects={len(np.unique(g))}")
head: X=(3321, 77)  raw=(3321, 100, 6)  dropped=0
chest: X=(3325, 77)  raw=(3325, 100, 6)  dropped=1
waist: X=(3324, 77)  raw=(3324, 100, 6)  dropped=2
wrist: X=(3320, 77)  raw=(3320, 100, 6)  dropped=6
thigh: X=(3319, 77)  raw=(3319, 100, 6)  dropped=7
ankle: X=(3319, 77)  raw=(3319, 100, 6)  dropped=6
cache written to data/cache.npz
head    X=(3321, 77)  fall_rate=0.553  subjects=17
chest   X=(3325, 77)  fall_rate=0.554  subjects=17
waist   X=(3324, 77)  fall_rate=0.554  subjects=17
wrist   X=(3320, 77)  fall_rate=0.554  subjects=17
thigh   X=(3319, 77)  fall_rate=0.554  subjects=17
ankle   X=(3319, 77)  fall_rate=0.554  subjects=17

4. Models, and a harness that does not leak¶

I lead with a Random Forest (200 trees) and a support vector machine (RBF kernel). Both are strong, well-understood classifiers for tabular features and both are defensible in a viva, which matters more here than chasing the last fraction of a percent. A small convolutional network appears later, but only as a deliberate counter-example, not as the method I am recommending.

The harness is where this assignment is won or lost. Leave-one-subject-out cross-validation holds out one entire volunteer per fold, trains on the other sixteen, and repeats seventeen times. The subtle mistake it is easy to make — and the exact mistake Assignment 1 criticised in the literature — is to scale or select features using the whole dataset before splitting. That leaks information about the held-out subject into training. I prevent it structurally: the scaler lives inside an sklearn Pipeline, so it is re-fit on each fold's training data only and then applied to the held-out subject. No fit touches the full feature matrix anywhere in the evaluation path. Two cells in Section 9 do fit on all data, for the feature-importance chart and the PCA projection; both are descriptive figures, neither contributes to any reported score, and I flag them here rather than leave the claim absolute.

Hyperparameters are fixed in advance, not selected from the data. The Random Forest uses 200 trees and the SVM an RBF kernel with C = 10, both chosen before seeing any result and never tuned by inspecting scores. Because no hyperparameter is selected from the data, there is no inner selection loop to leak through and no nested cross-validation is required. Had I tuned, the search would have to sit inside each training fold — the same principle the scaler already follows.

In [8]:
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import RandomForestClassifier
from sklearn.svm import SVC
from sklearn.model_selection import (LeaveOneGroupOut, StratifiedKFold,
                                     cross_val_predict)
from sklearn.metrics import (accuracy_score, precision_score, recall_score,
                             f1_score, confusion_matrix)

def make_model(kind):
    # the scaler is INSIDE the pipeline -> fit on train fold only, applied to held-out subject
    if kind == "rf":
        clf = RandomForestClassifier(n_estimators=200, random_state=0, n_jobs=-1)
    else:
        clf = SVC(kernel="rbf", C=10, gamma="scale")
    return Pipeline([("scale", StandardScaler()), ("clf", clf)])

def loso_predictions(X, y, g, kind):
    # out-of-fold prediction for every trial, each from a model blind to that subject
    return cross_val_predict(make_model(kind), X, y, groups=g,
                             cv=LeaveOneGroupOut(), n_jobs=-1)

def metric_row(y, pred):
    tn, fp, fn, tp = confusion_matrix(y, pred).ravel()
    return dict(accuracy=accuracy_score(y, pred),
                precision=precision_score(y, pred, zero_division=0),
                recall=recall_score(y, pred, zero_division=0),          # fall sensitivity
                specificity=tn / (tn + fp) if (tn + fp) else 0.0,
                f1=f1_score(y, pred, zero_division=0))

5. The cost of an honest test: k-fold versus LOSO¶

Before comparing placements, I show what the choice of validation is worth, on a single location, so the rest of the notebook can be read in that light. I run the wrist under both the subject-blind 10-fold cross-validation the earlier papers used and the subject-independent LOSO I argue for. No prior study on this dataset has put the two side by side.

In [9]:
loc = "wrist"
X, y, g = feats[loc]
loso_pred = loso_predictions(X, y, g, "rf")
kf_pred = cross_val_predict(make_model("rf"), X, y,
                            cv=StratifiedKFold(10, shuffle=True, random_state=0), n_jobs=-1)

comp = pd.DataFrame({"10-fold (subject-blind)": metric_row(y, kf_pred),
                     "LOSO (subject-independent)": metric_row(y, loso_pred)}).T
display(comp.round(4))
gap = comp.loc["10-fold (subject-blind)", "accuracy"] - comp.loc["LOSO (subject-independent)", "accuracy"]
print(f"accuracy inflation from subject leakage: {gap*100:.1f} points")
accuracy precision recall specificity f1
10-fold (subject-blind) 0.9780 0.9752 0.9853 0.9689 0.9803
LOSO (subject-independent) 0.9672 0.9616 0.9799 0.9514 0.9706
accuracy inflation from subject leakage: 1.1 points

The gap is real and in the expected direction: the subject-blind protocol reports the higher number, because it can recognise the individual rather than the fall. On this dataset the inflation is modest — around a point of accuracy — rather than the ten to fifteen points seen in some activity-recognition work. That is itself a finding worth stating plainly: these simulated falls are vigorous and separate cleanly even across people, so leakage flatters the result less than it might elsewhere. It does not flatter it by nothing, and the direction is the same one the literature warns about.

6. The placement comparison under LOSO¶

Now the main experiment: each of the six locations, evaluated under LOSO with both the Random Forest and the SVM, reporting the full set of metrics. Because falls are the high-cost class, recall on the fall class (sensitivity) and F1 matter more than raw accuracy. Two views of the same folds are computed and each has a job: pooled out-of-fold predictions, where every trial is scored by the model that did not see its subject, give the aggregate metrics and the confusion matrix; per-subject F1 — seventeen values per location — is what the paired significance tests in Section 7 consume. The pooled numbers are the headline; the per-subject spread is the evidence that the gaps are or are not real.

In [10]:
def per_subject_f1(X, y, g, kind):
    scores, subs = [], []
    for tr, te in LeaveOneGroupOut().split(X, y, g):
        m = make_model(kind).fit(X[tr], y[tr])
        scores.append(f1_score(y[te], m.predict(X[te]), zero_division=0))
        subs.append(np.unique(g[te])[0])
    return np.array(scores), np.array(subs)

rows, persubj = [], {}
t0 = time.time()
for kind in ["rf", "svm"]:
    for loc in LOC_ORDER:
        X, y, g = feats[loc]
        pred = loso_predictions(X, y, g, kind)
        m = metric_row(y, pred); m["model"] = kind.upper(); m["location"] = loc
        rows.append(m)
        if kind == "rf":
            s, _ = per_subject_f1(X, y, g, kind)
            persubj[loc] = s
print("LOSO sweep done in %.0fs" % (time.time() - t0))

results = pd.DataFrame(rows)[["model", "location", "accuracy", "precision",
                              "recall", "specificity", "f1"]]
results = results.sort_values(["model", "f1"], ascending=[True, False]).reset_index(drop=True)
display(results.round(4))
LOSO sweep done in 130s
model location accuracy precision recall specificity f1
0 RF ankle 0.9928 0.9929 0.9940 0.9912 0.9935
1 RF chest 0.9919 0.9913 0.9940 0.9892 0.9927
2 RF thigh 0.9916 0.9892 0.9956 0.9865 0.9924
3 RF waist 0.9904 0.9902 0.9924 0.9879 0.9913
4 RF head 0.9792 0.9722 0.9908 0.9649 0.9814
5 RF wrist 0.9672 0.9616 0.9799 0.9514 0.9706
6 SVM thigh 0.9952 0.9951 0.9962 0.9939 0.9956
7 SVM chest 0.9934 0.9924 0.9957 0.9906 0.9940
8 SVM ankle 0.9928 0.9929 0.9940 0.9912 0.9935
9 SVM waist 0.9922 0.9892 0.9967 0.9865 0.9930
10 SVM head 0.9840 0.9874 0.9837 0.9845 0.9856
11 SVM wrist 0.9768 0.9767 0.9815 0.9710 0.9791
In [11]:
# ranked summary, Random Forest, by fall-class F1
rf_res = results[results.model == "RF"].set_index("location").loc[LOC_ORDER]
best = rf_res.f1.idxmax()
print("best location (RF, fall F1): %s = %.3f" % (best, rf_res.f1.max()))
print("wrist (RF, fall F1):         %.3f" % rf_res.loc["wrist", "f1"])
for other in ["waist", "thigh"]:
    print("  wrist - %-5s = %+.3f F1" % (other, rf_res.loc["wrist", "f1"] - rf_res.loc[other, "f1"]))
best location (RF, fall F1): ankle = 0.993
wrist (RF, fall F1):         0.971
  wrist - waist = -0.021 F1
  wrist - thigh = -0.022 F1

7. Are the placement gaps real, or noise?¶

Seventeen per-subject F1 scores per location is exactly the paired sample a significance test wants. The scores are not normally distributed and the sample is small, so I use the Wilcoxon signed-rank test (paired, non-parametric) rather than a t-test, and I add a bootstrap 95% confidence interval on the mean paired difference. This answers the question the bar chart alone cannot: when the wrist trails the waist and thigh, is that a real effect or the luck of which subjects landed where?

In [12]:
def paired_ci(a, b, n_boot=10000, seed=0):
    rng = np.random.default_rng(seed)
    diff = a - b
    boot = np.array([rng.choice(diff, len(diff), replace=True).mean() for _ in range(n_boot)])
    return diff.mean(), np.percentile(boot, 2.5), np.percentile(boot, 97.5)

print("Wrist versus the two leading trunk/thigh locations (RF, per-subject fall F1):\n")
for other in ["waist", "thigh"]:
    a, b = persubj["wrist"], persubj[other]
    stat, p = wilcoxon(a, b)
    mean_d, lo, hi = paired_ci(a, b)
    print(f"wrist vs {other:5s}:  mean diff = {mean_d:+.3f} F1  "
          f"95% CI [{lo:+.3f}, {hi:+.3f}]  Wilcoxon p = {p:.4f}")
Wrist versus the two leading trunk/thigh locations (RF, per-subject fall F1):

wrist vs waist:  mean diff = -0.021 F1  95% CI [-0.029, -0.013]  Wilcoxon p = 0.0008
wrist vs thigh:  mean diff = -0.022 F1  95% CI [-0.029, -0.016]  Wilcoxon p = 0.0005

The gap between the wrist and the best placements is small in absolute terms — around two F1 points — but the confidence interval sits clear of zero and the Wilcoxon test rejects the null, so the wrist genuinely is the weaker location, not just unlucky.

It would be inconsistent to test the wrist gap and then merely assert that the top four locations are a tie — that eyeballing is exactly what Assignment 1 criticises. So I test it. Following Demšar's (2006) recommended procedure for comparing methods over a common sample, I use the Friedman test as an omnibus across all six locations, then Wilcoxon signed-rank for the pairs within the leading group, with a Holm correction for the six comparisons.

In [13]:
from itertools import combinations
from scipy.stats import friedmanchisquare

def holm(pvals):
    # Holm-Bonferroni step-down adjustment
    p = np.asarray(pvals, float); m = len(p); adj = np.empty(m); running = 0.0
    for i, idx in enumerate(np.argsort(p)):
        running = max(running, (m - i) * p[idx]); adj[idx] = min(running, 1.0)
    return adj

# omnibus: do the six locations differ at all?
stat, p_omni = friedmanchisquare(*[persubj[l] for l in LOC_ORDER])
print(f"Friedman across six locations: chi2 = {stat:.2f}, p = {p_omni:.2e}\n")

# leading group, pairwise, corrected for the six comparisons
LEAD = ["chest", "waist", "thigh", "ankle"]
rows = []
for a, b in combinations(LEAD, 2):
    _, p = wilcoxon(persubj[a], persubj[b])
    md_, lo, hi = paired_ci(persubj[a], persubj[b])
    rows.append(dict(pair=f"{a} vs {b}", mean_diff=md_, ci_low=lo, ci_high=hi, p_raw=p))
lead = pd.DataFrame(rows)
lead["p_holm"] = holm(lead.p_raw.values)
display(lead.round(4))
print("significant after Holm:", int((lead.p_holm < 0.05).sum()), "of", len(lead))
Friedman across six locations: chi2 = 36.43, p = 7.78e-07

pair mean_diff ci_low ci_high p_raw p_holm
0 chest vs waist 0.0012 -0.0028 0.0060 0.9164 1.0
1 chest vs thigh 0.0002 -0.0040 0.0050 0.8753 1.0
2 chest vs ankle -0.0005 -0.0041 0.0032 0.8753 1.0
3 waist vs thigh -0.0010 -0.0035 0.0018 0.3850 1.0
4 waist vs ankle -0.0017 -0.0063 0.0023 0.5748 1.0
5 thigh vs ankle -0.0007 -0.0054 0.0037 0.9594 1.0
significant after Holm: 0 of 6

The two results define the picture precisely. The Friedman test is decisive (p below 1e-6): the six locations are not all equivalent — but that omnibus is driven by the wrist and, to a lesser extent, the head. Within the leading group of chest, waist, thigh and ankle, none of the six pairwise differences survives the Holm correction; every adjusted p-value is 1.0 and every mean difference is under two thousandths of an F1 point. So the top four locations are statistically indistinguishable on this data, and the earlier remark that ankle "edged" the cluster is noise, not a ranking. The disadvantage of the wrist, by contrast, is real and it is small — both facts belong in the conclusion.

8. Why not deep learning here — a demonstration, not a strawman¶

Assignment 1 argued that a data-hungry network is the wrong tool for 17 subjects. It would be too easy to assert that and move on, so I show it. Below is a small 1-D convolutional network trained on the raw windowed signal of the best location, under the same LOSO protocol, with normalisation fit inside each fold. If the argument is right, it should not beat the Random Forest.

In [14]:
import torch, torch.nn as nn
torch.manual_seed(0)

class SmallCNN(nn.Module):
    def __init__(self, ch=6):
        super().__init__()
        self.net = nn.Sequential(
            nn.Conv1d(ch, 32, 5, padding=2), nn.ReLU(), nn.MaxPool1d(2),
            nn.Conv1d(32, 64, 5, padding=2), nn.ReLU(), nn.AdaptiveAvgPool1d(1),
            nn.Flatten(), nn.Dropout(0.3), nn.Linear(64, 2))
    def forward(self, x): return self.net(x)

def cnn_fold(Wtr, ytr, Wte, epochs=12):
    mu, sd = Wtr.mean((0, 1), keepdims=True), Wtr.std((0, 1), keepdims=True) + 1e-6
    Wtr, Wte = (Wtr - mu) / sd, (Wte - mu) / sd
    Xtr = torch.tensor(Wtr.transpose(0, 2, 1).astype("float32"))
    Xte = torch.tensor(Wte.transpose(0, 2, 1).astype("float32"))
    yt = torch.tensor(ytr).long()
    m = SmallCNN(); opt = torch.optim.Adam(m.parameters(), lr=1e-3); lf = nn.CrossEntropyLoss()
    m.train()
    for _ in range(epochs):
        for i in range(0, len(Xtr), 64):
            opt.zero_grad(); lf(m(Xtr[i:i+64]), yt[i:i+64]).backward(); opt.step()
    m.eval()
    with torch.no_grad():
        return m(Xte).argmax(1).numpy()

W, yW, gW = raws[best]
t0 = time.time()
cnn_pred, cnn_true = [], []
for tr, te in LeaveOneGroupOut().split(W, yW, gW):
    cnn_pred.append(cnn_fold(W[tr], yW[tr], W[te])); cnn_true.append(yW[te])
cnn_pred, cnn_true = np.concatenate(cnn_pred), np.concatenate(cnn_true)
cnn_m = metric_row(cnn_true, cnn_pred)
print("1-D CNN on %s (LOSO): acc=%.3f  fall F1=%.3f   [%.0fs]" % (
    best, cnn_m["accuracy"], cnn_m["f1"], time.time() - t0))
print("Random Forest on %s (LOSO): acc=%.3f  fall F1=%.3f" % (
    best, rf_res.loc[best, "accuracy"], rf_res.loc[best, "f1"]))
1-D CNN on ankle (LOSO): acc=0.906  fall F1=0.911   [46s]
Random Forest on ankle (LOSO): acc=0.993  fall F1=0.993

The network lands several F1 points below the Random Forest on the same location and the same folds. The comparison is symmetric in the way that matters: neither model received any hyperparameter search. The network trains for a fixed twelve epochs at a fixed learning rate, and the forest is equally untuned at 200 trees. The result therefore shows that with 17 subjects and this feature set the classical model reaches a level this network does not — it does not show that no network could, and I would not claim that. That is the demonstration behind Assignment 1's argument, offered as a measured result rather than a rhetorical one.

9. Figures¶

The first figure is the headline: fall-class F1 per location under LOSO, with error bars across the seventeen subjects, so the answer — which location, and how confidently — is legible at a glance. It is followed by the confusion matrix of the best location, the Random Forest's feature importances, and a PCA projection of the feature space coloured by class.

In [15]:
# Figure 1 — headline: per-location fall F1 (mean +/- SD across 17 subjects), RF
fig, ax = plt.subplots(figsize=(7.5, 4.2))
means = [persubj[l].mean() for l in LOC_ORDER]
sds = [persubj[l].std() for l in LOC_ORDER]
bars = ax.bar(LOC_ORDER, means, yerr=sds, capsize=4,
              color=[PALETTE[l] for l in LOC_ORDER], edgecolor="white")
ax.set_ylim(0.9, 1.0); ax.set_ylabel("fall-class F1 (LOSO)")
ax.set_title("Fall-detection F1 by sensor location, subject-independent")
for l, m in zip(LOC_ORDER, means):
    ax.text(l, m + 0.002, f"{m:.3f}", ha="center", va="bottom", fontsize=9)
ax.axhline(max(means), ls="--", lw=0.8, color="#888")
plt.tight_layout(); plt.savefig("output/fig_placement.png", dpi=130); plt.show()
No description has been provided for this image
In [16]:
# Figure 2 — confusion matrix of the best location (RF, LOSO out-of-fold predictions)
Xb, yb, gb = feats[best]
pb = loso_predictions(Xb, yb, gb, "rf")
cm = confusion_matrix(yb, pb)
fig, ax = plt.subplots(figsize=(4.2, 3.8))
im = ax.imshow(cm, cmap="Blues")
ax.set_xticks([0, 1], ["ADL", "fall"]); ax.set_yticks([0, 1], ["ADL", "fall"])
ax.set_xlabel("predicted"); ax.set_ylabel("true")
ax.set_title(f"Confusion matrix — {best} (RF, LOSO)")
for i in range(2):
    for j in range(2):
        ax.text(j, i, cm[i, j], ha="center", va="center",
                color="white" if cm[i, j] > cm.max() / 2 else "black", fontsize=12)
ax.grid(False); plt.tight_layout(); plt.savefig("output/fig_confusion.png", dpi=130); plt.show()
No description has been provided for this image
In [17]:
# Figure 3 — Random Forest feature importances on the best location (top 15)
rf_full = make_model("rf").fit(Xb, yb)      # descriptive only — fitted on all data, contributes to no reported metric
imp = pd.Series(rf_full.named_steps["clf"].feature_importances_, index=feat_names)
top = imp.sort_values(ascending=False).head(15)[::-1]
fig, ax = plt.subplots(figsize=(6.8, 4.4))
ax.barh(top.index, top.values, color="#2bb2a8")
ax.set_title(f"Most informative features — {best}")
ax.set_xlabel("Random Forest importance")
plt.tight_layout(); plt.savefig("output/fig_importance.png", dpi=130); plt.show()
No description has been provided for this image
In [18]:
# Figure 4 — PCA of the best location's feature space, coloured by class
# descriptive only — fitted on all data, contributes to no reported metric
from sklearn.decomposition import PCA
Z = PCA(n_components=2, random_state=0).fit_transform(StandardScaler().fit_transform(Xb))
fig, ax = plt.subplots(figsize=(5.6, 4.6))
for lab, name, c in [(0, "ADL", "#4c9be8"), (1, "fall", "#e8734c")]:
    ax.scatter(Z[yb == lab, 0], Z[yb == lab, 1], s=6, alpha=0.4, label=name, color=c)
# a few real trials sit far out (violent falls) and would squash the cluster into a
# corner; clip the view to the dense region and report how many points fall outside it.
xlo, xhi = np.percentile(Z[:, 0], 1) - 2, np.percentile(Z[:, 0], 99) + 2
ylo, yhi = np.percentile(Z[:, 1], 1) - 2, np.percentile(Z[:, 1], 99) + 2
hidden = int(((Z[:, 0] < xlo) | (Z[:, 0] > xhi) | (Z[:, 1] < ylo) | (Z[:, 1] > yhi)).sum())
ax.set_xlim(xlo, xhi); ax.set_ylim(ylo, yhi)
ax.set_xlabel("PC1"); ax.set_ylabel("PC2"); ax.legend()
ax.set_title(f"Feature space, {best} (PCA, coloured by class)")
ax.text(0.99, 0.01, f"{hidden} real outlier trials outside axes", transform=ax.transAxes,
        ha="right", va="bottom", fontsize=8, color="#666")
plt.tight_layout(); plt.savefig("output/fig_pca.png", dpi=130); plt.show()
print(f"axes clipped to [{xlo:.1f},{xhi:.1f}] x [{ylo:.1f},{yhi:.1f}]; {hidden} real trials off-plot")
No description has been provided for this image
axes clipped to [-7.8,10.6] x [-8.7,5.5]; 9 real trials off-plot

10. Findings and conclusions¶

Read under honest validation, the placement picture is clear and a little more interesting than the published rankings suggest.

  • Several trunk and lower-body locations tie for the strongest placement. The chest, waist, thigh and ankle cluster within two thousandths of an F1 point of each other under LOSO, and a Friedman omnibus followed by Holm-corrected Wilcoxon tests (Section 7) confirms the four are statistically indistinguishable — the ankle's nominal lead is noise, not a ranking — while the head sits a step behind. This is consistent with Özdemir (2016) and Ntanasis et al. (2017), who named the waist the best placement; the waist is squarely in the leading group, now established on the 17-subject release and without subject leakage, which none of them did.
  • The wrist is the weakest location, again as the literature found. The new part is the size of the gap: about two F1 points behind the best placements. The Wilcoxon test confirms the gap is real, and the confidence interval confirms it is small.
  • Subject leakage inflates results on this dataset, but only modestly — around a point of accuracy on the wrist between 10-fold and LOSO. The vigorous, well-separated simulated falls are the reason; the inflation is smaller than in general activity recognition, but it is present and in the direction the literature warns of.
  • The deep network underperforms the Random Forest on the same location and folds, which is the demonstration behind Assignment 1's methodological argument rather than an assumption.

The decision this notebook was built to make: because the wrist sits within about two points of the best trunk or thigh placement under subject-independent testing, the everyday acceptability of a wristband is a defensible trade for a small, real accuracy cost — the best placement for a product need not be the best placement on paper. The argument is not that the wrist is as good; the significance test says it is not. It is that the honest gap is small enough for comfort and adherence to legitimately outweigh it, which is precisely the judgement the earlier k-fold numbers were too inflated to support.

Further work. The natural next step is to test whether a lightweight per-user calibration — a few of the wearer's own ADLs folded into training — closes the wrist's remaining two points, since the leakage analysis suggests a good part of the gap is person-specific. A second thread is cost-sensitive training, given that a missed fall is far more expensive than a false alarm and the current models optimise neither for that asymmetry.