Reproducing Masked Global Weighted Pearson with Leakage-Safe Validation

Hello! Last lesson established how to read one ordered sequence at a time without loading the full Parquet dataset. That gives us a bounded data pipeline, but it does not yet tell us whether a model is improving. For that, we need two things that agree with the competition: a correct scorer and a validation split that does not reveal the future.

This lesson takes about 40–45 minutes. You will derive a masked weighted Pearson score, implement a streaming version that works with sequence-sized batches, and design chronological validation. One important boundary: the supplied materials explain weighted correlation and time-series splitting, but do not contain the competition’s official scoring code. We can implement and test the stated metric structure here; exact reproduction must ultimately be checked against the organizer’s specification, particularly how its two outputs are combined.

What “masked,” “weighted,” and “global” each change

Pearson correlation measures whether predictions and targets vary together after each has been centered on its mean. Weighted Pearson replaces ordinary means with weighted means. The NIST reference gives the weighted means, covariance, and correlation in one place.

Suppose an eligible scored entry has target $y_i$, prediction ${\hat{y}}_i$, and nonnegative weight $w_i$. Define

The weighted Pearson correlation is

A mask decides which entries are eligible before calculating either mean. An unscored time step must not affect the means, variances, or numerator. A zero-weight entry likewise contributes nothing. Replacing masked targets with zero and then calculating correlation over the full array is wrong: those zeros change the centering.

“Global” describes where the calculation takes place. In a pooled global score, collect every eligible entry across validation sequences, calculate one pair of weighted means, then calculate one correlation. Computing a correlation separately for each sequence and averaging the results answers a different question. Likewise, averaging fold scores is a useful development summary, but it is not the same calculation as scoring the pooled predictions from those folds.

There is an additional ambiguity with two predictions per step. “Global” might mean pooling both output channels into one correlation, or calculating one global correlation per channel and combining the two scores. Those choices can produce different rankings. Before treating any local number as the competition score, confirm from its scoring contract:

• which time steps and output channels are eligible;
• whether weights apply per step, per channel, or another way;
• whether the two channels are pooled or scored separately, and any channel-combination rule;
• what happens when a scored target or prediction is non-finite, or a variance is zero.

The implementation below deliberately adopts one explicit, provisional contract: mask and weights have the same $(T,2)$ shape as targets and predictions, and all eligible entries from both channels are pooled. Do not silently substitute it for a different official rule.

A scorer that stays streaming

The previous lesson’s loader yielded complete sequences. To score them without concatenating every validation tensor, we need quantities that can be accumulated across sequences. Expanding the centered sums gives

with analogous expressions

After each sequence, we retain only these running sums, not its full predictions. The following NumPy implementation enforces the provisional contract. It accepts masked-out target values that are non-finite, but rejects non-finite values where they would be scored.

import numpy as np


class PooledWeightedPearson:
    """One correlation across all eligible entries of (T, 2) arrays."""

    def __init__(self):
        self.W = 0.0
        self.wy = 0.0
        self.wp = 0.0
        self.wyy = 0.0
        self.wpp = 0.0
        self.wyp = 0.0

    def update(self, target, prediction, mask, weight):
        y = np.asarray(target, dtype=np.float64)
        p = np.asarray(prediction, dtype=np.float64)
        m = np.asarray(mask)
        w = np.asarray(weight, dtype=np.float64)

        if y.ndim != 2 or y.shape[1] != 2:
            raise ValueError("Expected targets with shape (T, 2)")
        if not (p.shape == m.shape == w.shape == y.shape):
            raise ValueError("All four arrays must have identical shapes")
        if m.dtype != np.bool_:
            raise ValueError("Mask must contain booleans")
        if not np.isfinite(w).all() or (w < 0).any():
            raise ValueError("Weights must be finite and nonnegative")

        use = m & (w > 0)
        y, p, w = y[use], p[use], w[use]
        if not np.isfinite(y).all() or not np.isfinite(p).all():
            raise ValueError("Non-finite value in a scored entry")

        self.W += np.sum(w)
        self.wy += np.sum(w * y)
        self.wp += np.sum(w * p)
        self.wyy += np.sum(w * y * y)
        self.wpp += np.sum(w * p * p)
        self.wyp += np.sum(w * y * p)

    def result(self):
        if self.W <= 0:
            raise ValueError("No positive-weight scored entries")

        covariance = self.wyp - self.wy * self.wp / self.W
        target_spread = self.wyy - self.wy * self.wy / self.W
        prediction_spread = self.wpp - self.wp * self.wp / self.W
        if target_spread <= 0 or prediction_spread <= 0:
            raise ValueError("Correlation undefined for zero spread")

        return covariance / np.sqrt(target_spread * prediction_spread)

For a loader that also yields mask and weight, call update once per validation sequence. The earlier demonstration loader yielded only x and y; it would need to read the actual mask and weight fields, or derive them exactly as the scoring rules specify. Do not invent an all-true mask or unit weights merely to make the code run.

A compact check exposes two common mistakes:

target = np.array([[0., 0.], [1., 1.], [2., 2.]])
prediction = np.array([[0., 0.], [1., 99.], [2., 2.]])
mask = np.array([[True, True], [True, False], [True, True]])
weight = np.ones((3, 2))

score = PooledWeightedPearson()
score.update(target, prediction, mask, weight)
assert np.isclose(score.result(), 1.0)

The conspicuously bad prediction is masked, so it cannot change the score. Unmask it, and the result should change. Also check that calling update on two parts of an array gives approximately the same result as calling it once on the whole array. That tests the intended global accumulation, not merely the arithmetic within one sequence.

This raw-sum formulation uses float64, but subtracting large, nearly equal sums can still lose precision when values have huge offsets and tiny variation. If that occurs in the real data, use a numerically stable weighted-moment accumulator or a two-pass computation, and compare both against a small, directly centered reference. Neither precision fix changes the definition of the score.

Put validation later in time

A metric implementation can be mathematically perfect and still report a misleading result if its validation data leaks into training. Watch the short demonstration of a time-series split.

An expanding-window procedure trains a fresh model on the available earlier periods for each fold and evaluates it on the next held-out period. “Fresh” matters: starting fold two with weights already trained on fold one’s later validation period can invalidate fold two’s training boundary. A sliding window instead limits how far back the training history extends; either can respect chronology.

The scikit-learn documentation provides a precise implementation of expanding training indices, including an optional gap.

For this competition, create split indices over ordered time blocks, not arbitrary rows returned by the Parquet loader. For example, if each complete sequence belongs to exactly one block, a planning sketch is:

import numpy as np
from sklearn.model_selection import TimeSeriesSplit

periods = np.array(sorted(available_training_periods))
splitter = TimeSeriesSplit(n_splits=3, test_size=2, gap=1)

for train_idx, val_idx in splitter.split(periods):
    train_periods = set(periods[train_idx])
    val_periods = set(periods[val_idx])

    assert max(train_periods) < min(val_periods)

    # Select complete sequences belonging to each set of periods.
    # Fit preprocessing and a NEW model on train sequences only.
    # Predict val sequences in their permitted chronological order.
    # Feed predictions, targets, masks, and weights to a NEW scorer.

Here gap=1 excludes one indexed period between each training and validation window. It does not automatically solve every leakage problem. If a target attached to a training feature time depends on observations extending into the validation window, purge training examples whose label-information interval crosses the boundary. Choose the gap using the actual target horizon and data semantics, not a conventional number.

Similarly, TimeSeriesSplit knows nothing about seq_id. If one complete sequence spans a proposed boundary, do not assign its early rows to training and late rows to validation while claiming a sequence-level holdout. Move the boundary or exclude the crossing sequence, according to the intended deployment test. Check identifiers and time ranges after selecting row groups, not just the printed period indices.

Make the score reflect the deployed prediction path

A safe split is more than an ordering assertion. In each fold, fit feature normalization, imputation statistics, target transforms, and model parameters only on the training portion. Apply the fitted transforms unchanged to validation. Features computed using later rows—such as a centered rolling window—can leak even when the split itself is chronological.

For a stateful sequence model, validation must also follow its inference contract. Reset state at each true independent sequence boundary. At a validation step, provide only features that would be available at that step; do not feed validation targets back as inputs unless the competition explicitly permits them. If the rules permit a feature-only warm-up prefix, it may update hidden state without contributing to the metric; the evaluation mask still controls which later predictions are scored. These choices should be checked against the same competition contract as the metric.

Maintain two distinct evaluation roles:

1. Use chronological development folds to compare model and preprocessing choices. Report each fold’s masked weighted Pearson score so a weak time period is visible, rather than hiding it in a mean.
2. Reserve a final, later holdout for a single assessment after those choices are settled. Do not repeatedly tune on that holdout.

Within one validation fold, accumulate the scorer globally over all its eligible sequence entries. If the competition defines a single pooled score across an evaluation set, do not replace it with a mean of per-sequence correlations. Conversely, a mean of fold scores may summarize stability, but label it as such—not as an exact reproduction of the competition calculation.

Finally, perform a parity check before relying on your local metric: run the same small synthetic predictions, masks, targets, and weights through both your implementation and the official scorer, if one is available. Include cases with masked extreme predictions, unequal weights, and two channels with different behavior. If they disagree, resolve eligibility, channel aggregation, or zero-variance rules before training a larger model. Without that check—or an unambiguous official formula—the implementation above remains a well-defined candidate, not a verified replica.

Takeaways

Masked weighted Pearson first selects eligible entries, then calculates weighted means and centered covariance. A global score pools those entries at the level specified by the scoring contract; streaming sums let you do that without keeping all validation predictions in memory. Chronological folds, intact sequence boundaries, training-only preprocessing, and a faithful stateful inference path make the resulting number meaningful. The outstanding verification step is to match the competition’s exact two-output aggregation and edge-case rules.

Next, you will implement a GRU that maintains state within a sequence, resets it between independent sequences, and emits both required predictions at each step—the model that this validation pipeline will evaluate.