mkdir -p ~/.kaggle
python3 -c "import json;print(json.load(open('kaggle.json'))['key'],end='')" > ~/.kaggle/access_token
chmod 600 ~/.kaggle/access_token

%matplotlib inline
import time, warnings; warnings.filterwarnings('ignore')   # keep output clean
import numpy as np, pandas as pd                            # numerics + dataframes
import matplotlib.pyplot as plt                             # static plots (embed in HTML+PDF)
plt.rcParams['figure.dpi'] = 120                            # crisp figures
RANDOM_STATE = 0                                            # single seed used everywhere
np.random.seed(RANDOM_STATE)                                # reproducible sampling
NEG_WORD, POS_WORD = 'benign', 'attack'                      # class names (overridden by some loaders)

import os, glob, re
# RBA (Wiefling, Jørgensen, Thunem & Lo Iacono, 2022): 31.3M login attempts to a large online
# service. The public Kaggle release is a SYNTHESIZED version of the original production traffic
# (re-generated so customers cannot be re-identified), so it mirrors real behaviour but is not
# itself the raw log,
# labelled `Is Attack IP` (login from a known-attacker IP) and `Is Account Takeover`. A risk-based-
# authentication / credential-attack detection problem. Self-contained Kaggle download (~9 GB).
os.environ.setdefault('KAGGLE_KEY', open(os.path.expanduser('~/.kaggle/access_token')).read().strip())
RBA_DIR = '/tmp/kg_rba'; os.makedirs(RBA_DIR, exist_ok=True)
hits = glob.glob(RBA_DIR + '/**/rba-dataset.csv', recursive=True)
if not hits:
    import kaggle; kaggle.api.authenticate()
    print('downloading RBA (~9 GB, one-time)...')
    kaggle.api.dataset_download_files('dasgroup/rba-dataset', path=RBA_DIR, unzip=True, quiet=True)
    hits = glob.glob(RBA_DIR + '/**/rba-dataset.csv', recursive=True)
f = hits[0]; assert os.path.exists(f), f'RBA csv not found: {f}'
# Read in CHUNKS (the User-Agent field contains commas, so proper CSV quoting matters, and the file
# is 9 GB). We DROP the IP Address, User ID and timestamp: `Is Attack IP` is literally a property of
# the IP, so those columns would be direct label leakage. We keep geo / ASN / client-string features
# and let the audit reveal how much the label rides on IP-derived metadata. Attack IPs are ~10% and
# time-clustered, so we KEEP EVERY attack row and SUBSAMPLE benign (bounded-sample rate, not base rate).
USE = ['Round-Trip Time [ms]','Country','Region','City','ASN','Browser Name and Version',
       'OS Name and Version','Device Type','Is Attack IP','Is Account Takeover']
frames = []
for ch in pd.read_csv(f, usecols=USE, chunksize=3_000_000, low_memory=False):
    ch.columns = [c.strip() for c in ch.columns]
    ap = ch['Is Attack IP'].astype(str).str.lower().isin(['true', '1'])
    neg = ch[~ap].sample(frac=0.12, random_state=0)
    frames.append(pd.concat([ch[ap], neg]))
df = pd.concat(frames, ignore_index=True)
assert len(df) >= 1_000_000, f'floor not met: {len(df):,}'
df['y'] = df['Is Attack IP'].astype(str).str.lower().isin(['true', '1']).astype(int)
# family = per-country grouping, so per-family recall exposes GEOGRAPHIC detection bias.
df['family'] = df['Country'].astype(str).fillna('??')
DROP = ['y', 'family', 'Is Attack IP', 'Is Account Takeover']
feat = [c for c in df.columns if c not in DROP]
from sklearn.preprocessing import LabelEncoder
X = df[feat].copy()
for c in X.select_dtypes(include='object').columns:
    X[c] = LabelEncoder().fit_transform(X[c].astype(str))
X = X.apply(pd.to_numeric, errors='coerce').replace([np.inf, -np.inf], np.nan).fillna(0.0)
X = X.loc[:, X.nunique() > 1]
X.columns = [re.sub(r'[^0-9A-Za-z_]+', '_', str(c)).strip('_') or f'f{i}'
             for i, c in enumerate(X.columns)]                 # LightGBM-safe names ([ms], spaces)
feat = list(X.columns)
y = df['y'].to_numpy(); family = df['family'].to_numpy()
print(f'loaded {len(df):,} login attempts x {len(feat)} features; '
      f'attack-IP rate (bounded sample) {y.mean():.4f}')

loaded 6,477,653 login attempts x 8 features; attack-IP rate (bounded sample) 0.4781

# --- EDA 1: class balance and the attack-family mix ---
fig, ax = plt.subplots(1, 2, figsize=(11, 4))
df['y'].map({0:NEG_WORD,1:POS_WORD}).value_counts().plot.bar(               # counts per class
    ax=ax[0], color=['#2a9d8f','#e76f51']); ax[0].set_yscale('log')
ax[0].set_title(f'Class balance ({NEG_WORD} vs {POS_WORD})'); ax[0].set_ylabel('records (log)')
df.loc[df.y==1,'family'].value_counts().head(8).plot.barh(                  # top attack families
    ax=ax[1], color='#e76f51'); ax[1].invert_yaxis(); ax[1].set_title('Top attack families')
plt.tight_layout(); plt.show()

# --- EDA 2: feature correlation + a 2-D PCA projection ---
from sklearn.preprocessing import StandardScaler                 # scale before PCA
from sklearn.decomposition import PCA
fig, ax = plt.subplots(1, 2, figsize=(12, 5))
topv = X[feat].var().sort_values().tail(12).index                # 12 highest-variance features
im = ax[0].imshow(X[topv].corr(), cmap='coolwarm', vmin=-1, vmax=1)  # correlation heatmap
ax[0].set_xticks(range(len(topv))); ax[0].set_xticklabels(topv, rotation=90, fontsize=7)
ax[0].set_yticks(range(len(topv))); ax[0].set_yticklabels(topv, fontsize=7)
ax[0].set_title('Feature correlation (top-variance)'); fig.colorbar(im, ax=ax[0], shrink=0.7)
samp = X.sample(min(5000, len(X)), random_state=RANDOM_STATE)     # subsample for a fast PCA
pc = PCA(n_components=2).fit_transform(StandardScaler().fit_transform(samp))
ys = y[samp.index]                                               # aligned labels for coloring
for lab,c in [(0,'#2a9d8f'),(1,'#e76f51')]:
    ax[1].scatter(pc[ys==lab,0], pc[ys==lab,1], s=4, alpha=0.4, color=c,
                  label={0:NEG_WORD,1:POS_WORD}[lab])
ax[1].set_title('PCA projection (2 components)'); ax[1].legend(); ax[1].set_xlabel('PC1'); ax[1].set_ylabel('PC2')
plt.tight_layout(); plt.show()

# --- Model comparison: four learners on the same held-out split ---
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import accuracy_score, roc_auc_score
import xgboost as xgb, lightgbm as lgb

# Stratified split keeps the class ratio in both halves.
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.25, random_state=RANDOM_STATE, stratify=y)
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
N_MATERIALIZED = len(y)                                          # the full corpus we loaded (see printed count)
# HONEST DISCLOSURE: we do NOT train on all N. We fit on a STRATIFIED subsample (<=120k) because
# these learners saturate long before then on this data. Every headline below is a SUBSAMPLE
# number, not a full-corpus number — saying otherwise would be the fabrication this course forbids.
if len(Xtr) > 120_000:
    Xtr, _, ytr, _ = train_test_split(Xtr, ytr, train_size=120_000, random_state=RANDOM_STATE,
                                      stratify=ytr)               # genuinely stratified, not random
MAJORITY_BASELINE = max(np.mean(yte), 1 - np.mean(yte))          # accuracy of 'always predict majority'
print(f'materialized {N_MATERIALIZED:,} rows | trained on {len(Xtr):,} (stratified subsample) | '
      f'held-out {len(yte):,}')
print(f'MAJORITY-CLASS BASELINE accuracy = {MAJORITY_BASELINE:.4f}  '
      f'(any model must beat THIS, not 0.5, to be interesting)')

models = {                                                        # four standard, diverse learners
    'LogisticRegression': make_pipeline(StandardScaler(), LogisticRegression(max_iter=300)),  # scaled!
    'RandomForest': RandomForestClassifier(n_estimators=60, n_jobs=-1, random_state=RANDOM_STATE),
    'XGBoost': xgb.XGBClassifier(n_estimators=80, max_depth=6, tree_method='hist', n_jobs=-1,
                                 eval_metric='logloss', random_state=RANDOM_STATE),
    'LightGBM': lgb.LGBMClassifier(n_estimators=80, n_jobs=-1, verbose=-1, random_state=RANDOM_STATE),
}
rows, fitted = [], {}
for name, m in models.items():                                    # fit + score each model
    t = time.perf_counter(); m.fit(Xtr, ytr); fitted[name] = m
    p = m.predict_proba(Xte)[:, 1]                                # positive-class probability on held-out
    rows.append({'model': name, 'accuracy': round(accuracy_score(yte, (p>0.5).astype(int)), 6),
                 'roc_auc': round(roc_auc_score(yte, p), 6),      # 6 dp: a 1.000000 is a red flag, not a win
                 'train_s': round(time.perf_counter()-t, 1)})
rows.append({'model': 'MajorityBaseline', 'accuracy': round(MAJORITY_BASELINE, 4),
             'roc_auc': 0.5, 'train_s': 0.0})            # show the baseline IN the ranking table
comparison = pd.DataFrame(rows).sort_values('roc_auc', ascending=False).reset_index(drop=True)
_ranked = comparison[comparison.model != 'MajorityBaseline']
best_name = _ranked.iloc[0]['model']; best = fitted[best_name]  # winner by ROC-AUC (excl. baseline)
print('best model:', best_name); comparison

materialized 6,477,653 rows | trained on 120,000 (stratified subsample) | held-out 1,619,414
MAJORITY-CLASS BASELINE accuracy = 0.5219  (any model must beat THIS, not 0.5, to be interesting)

best model: XGBoost

# --- Results for the best model: confusion, ROC, PR, importances, per-family recall ---
from sklearn.metrics import confusion_matrix, roc_curve, precision_recall_curve, recall_score
pb = best.predict_proba(Xte)[:, 1]; pred = (pb > 0.5).astype(int)
fig, ax = plt.subplots(1, 3, figsize=(15, 4))
# (1) confusion matrix
cm = confusion_matrix(yte, pred); ax[0].imshow(cm, cmap='Blues')
ax[0].set_title(f'{best_name}: confusion'); ax[0].set_xticks([0,1]); ax[0].set_yticks([0,1])
ax[0].set_xticklabels([NEG_WORD,POS_WORD]); ax[0].set_yticklabels([NEG_WORD,POS_WORD])
for (i,j),v in np.ndenumerate(cm): ax[0].text(j,i,f'{v:,}',ha='center',va='center')
# (2) ROC and PR curves
fpr,tpr,_ = roc_curve(yte, pb); prec,rec,_ = precision_recall_curve(yte, pb)
ax[1].plot(fpr,tpr,color='#264653'); ax[1].plot([0,1],[0,1],'--',c='grey')
ax[1].set_title(f'ROC (AUC={roc_auc_score(yte,pb):.4f})'); ax[1].set_xlabel('FPR'); ax[1].set_ylabel('TPR')
ax[2].plot(rec,prec,color='#e76f51'); ax[2].set_title('Precision-Recall'); ax[2].set_xlabel('recall'); ax[2].set_ylabel('precision')
plt.tight_layout(); plt.show()

# (3) feature importances + (4) per-attack-family recall
fig, ax = plt.subplots(1, 2, figsize=(13, 5))
imp, names = None, feat                                           # importances, robust to the scaled-LR pipeline
if hasattr(best, 'feature_importances_'):                          # tree models
    imp = best.feature_importances_; names = list(getattr(best, 'feature_names_in_', feat))[:len(imp)]
elif hasattr(best, 'named_steps') and 'logisticregression' in getattr(best, 'named_steps', {}):
    imp = np.abs(best.named_steps['logisticregression'].coef_[0]); names = feat  # LR pipeline
elif hasattr(best, 'coef_'):
    imp = np.abs(best.coef_[0]); names = feat
if imp is not None:
    pd.Series(imp, index=names[:len(imp)]).sort_values().tail(12).plot.barh(ax=ax[0], color='#264653')
ax[0].set_title(f'{best_name}: top importances / |coef|')
# Per-family recall, WORST-first so rare, hard classes are visible, not just the dominant floods.
fam_te = df.loc[Xte.index, 'family']
fr = {}
for fam, cnt in fam_te[yte==1].value_counts().items():
    if cnt < 5: continue                                          # need a few positives for a meaningful recall
    mask = (fam_te==fam).to_numpy(); fr[fam] = recall_score(yte[mask], pred[mask], zero_division=0)
srt = pd.Series(fr).sort_values()
show = pd.concat([srt.head(9), srt.tail(3)]) if len(srt) > 12 else srt   # worst 9 + best 3
show = show[~show.index.duplicated()]
show.plot.barh(ax=ax[1], color=['#e76f51' if v < 0.5 else '#2a9d8f' for v in show]); ax[1].set_xlim(0,1)
ax[1].set_title('Per-family recall (worst first; red < 0.5)')
plt.tight_layout(); plt.show()
# Operational numbers, not just figures: false-positive rate and the worst per-family recalls.
tn, fp = int(cm[0,0]), int(cm[0,1])
fpr_op = fp/(fp+tn) if (fp+tn) > 0 else float('nan')             # benign wrongly flagged @0.5
print(f'operational FALSE-POSITIVE RATE @0.5 = {fpr_op:.4f}  ({fp:,} benign flagged of {fp+tn:,})')
print('worst per-family recalls:', {k: round(v, 3) for k, v in srt.head(6).items()})

operational FALSE-POSITIVE RATE @0.5 = 0.2111  (178,440 benign flagged of 845,169)
worst per-family recalls: {'SM': 0.0, 'KH': 0.0, 'HN': 0.0, 'GR': 0.0, 'JM': 0.0, 'JO': 0.0}

# --- Validity audit: is the score real detection, or a data shortcut? ---
from sklearn.metrics import roc_auc_score
samp = X.sample(min(60_000, len(X)), random_state=1); ysamp = y[samp.index]
aucs = {}
for c in feat:                                                    # AUC of EACH feature alone
    col = samp[c].to_numpy(float)
    if col.std()==0: continue
    a = roc_auc_score(ysamp, col); aucs[c] = max(a, 1-a)          # direction-agnostic
best_auc = max(aucs.values()); best_col = max(aucs, key=aucs.get)
dup_rate = 1 - X.drop_duplicates().shape[0]/len(X)               # exact-duplicate feature rows (whole set)
# The statistic that actually inflates a held-out score is TRAIN/TEST CONTAMINATION: how many test
# rows are exact duplicates of a training row. Measure it directly on the split used above.
_trkeys = set(map(tuple, np.round(Xtr.to_numpy(), 6)))
_te = np.round(Xte.to_numpy(), 6)[:50_000]
contam = float(np.mean([tuple(r) in _trkeys for r in _te]))      # fraction of test rows seen in train
# Trust grade reflects BOTH failure modes and takes the WORSE of the two: a near-perfect single
# feature (shortcut) OR heavy train/test contamination each independently invalidate the headline.
_ga = 'F' if best_auc>=0.999 else 'D' if best_auc>=0.99 else 'C' if best_auc>=0.95 else 'B' if best_auc>=0.85 else 'A'
_gc = 'F' if contam>=0.5 else 'D' if contam>=0.3 else 'C' if contam>=0.15 else 'B' if contam>=0.05 else 'A'
grade = max(_ga, _gc)                                            # 'max' letter = worse grade (A best, F worst)
print(f'best single-feature AUC = {best_auc:.4f}  (feature: {best_col})')
print(f'   note: a near-1.0 single-feature AUC means this feature is *near-sufficient* (a shortcut),\n'
      f'   which may be legitimate signal OR an artifact — it is NOT the same as target leakage.')
print(f'exact-duplicate row rate (whole corpus) = {dup_rate:.3f}')
print(f'TRAIN/TEST exact-row contamination       = {contam:.3f}  (single-feat grade {_ga}, contam grade {_gc})')
print(f'==> data trust grade: {grade}   (worse of the two; F = shortcut and/or heavy contamination)')
s = pd.Series(aucs).sort_values().tail(15)
fig, ax = plt.subplots(figsize=(8,5))
s.plot.barh(ax=ax, color=['#e76f51' if v>=0.99 else '#457b9d' for v in s]); ax.axvline(0.5,ls='--',c='grey')
ax.set_xlim(0.5,1.0); ax.set_title('Single-feature ROC-AUC (red = near-perfect shortcut)'); ax.set_xlabel('AUC alone')
plt.tight_layout(); plt.show()

best single-feature AUC = 0.8095  (feature: ASN)
   note: a near-1.0 single-feature AUC means this feature is *near-sufficient* (a shortcut),
   which may be legitimate signal OR an artifact — it is NOT the same as target leakage.
exact-duplicate row rate (whole corpus) = 0.897
TRAIN/TEST exact-row contamination       = 0.737  (single-feat grade A, contam grade F)
==> data trust grade: F   (worse of the two; F = shortcut and/or heavy contamination)

# --- Ablation: SHOW the inflation empirically, don't just narrate it ---
from sklearn.base import clone
def _retrain_auc(Xa, ya):                                        # re-split, stratified-subsample, refit best family
    xtr, xte, ytr2, yte2 = train_test_split(Xa, ya, test_size=0.25, random_state=RANDOM_STATE, stratify=ya)
    if len(xtr) > 120_000:
        xtr, _, ytr2, _ = train_test_split(xtr, ytr2, train_size=120_000, random_state=RANDOM_STATE, stratify=ytr2)
    m = clone(best); m.fit(xtr, ytr2)
    return roc_auc_score(yte2, m.predict_proba(xte)[:, 1])
base_auc = roc_auc_score(yte, best.predict_proba(Xte)[:, 1])     # (0) the headline held-out AUC
Xdd = X.drop_duplicates(); ydd = y[Xdd.index]                    # (1) de-duplicated corpus
auc_dedup = _retrain_auc(Xdd, ydd)
auc_noshort = _retrain_auc(X.drop(columns=[best_col]), y) if best_col in X.columns else base_auc  # (2) drop shortcut
ablation = pd.DataFrame([
    {'setting': 'headline (as-is)',              'held_out_auc': round(base_auc, 6)},
    {'setting': f'de-duplicated ({1-len(Xdd)/len(X):.0%} rows removed)', 'held_out_auc': round(auc_dedup, 6)},
    {'setting': f'shortcut feature dropped ({best_col})', 'held_out_auc': round(auc_noshort, 6)},
])
print('Ablation — how much of the headline survives once each artifact is removed:')
ablation

Ablation — how much of the headline survives once each artifact is removed:

# --- Reproducibility & robustness ---
import sklearn
from sklearn.model_selection import StratifiedKFold, cross_val_score
print(f'seed={RANDOM_STATE} | numpy {np.__version__} | sklearn {sklearn.__version__} | '
      f'xgboost {xgb.__version__} | lightgbm {lgb.__version__}')
# 3-fold cross-validated ROC-AUC of the winning model (fresh clone, bounded subsample) -> mean +/- std.
from sklearn.base import clone
cvX, cvy = Xtr.iloc[:40_000], ytr[:40_000]
def _auc_scorer(est, Xv, yv):                                   # robust to xgboost's 2-col predict_proba
    p = est.predict_proba(Xv)
    p = p[:, 1] if getattr(p, 'ndim', 1) == 2 else p
    return roc_auc_score(yv, p)
try:
    cv = cross_val_score(clone(best), cvX, cvy,
                         cv=StratifiedKFold(3, shuffle=True, random_state=RANDOM_STATE),
                         scoring=_auc_scorer, error_score='raise')
    assert np.all(np.isfinite(cv)), 'non-finite CV folds'   # FAIL CLOSED: never narrate a NaN as evidence
    print(f'{best_name} 3-fold CV ROC-AUC = {cv.mean():.4f} +/- {cv.std():.4f}  '
          f'(mean +/- std across 3 stratified folds; a small std means a stable estimate on this split)')
except Exception as e:
    print(f'CV UNAVAILABLE ({type(e).__name__}: {str(e)[:60]}); rely on the single held-out AUC above — '
          f'we do NOT report a CV number we could not compute')

seed=0 | numpy 2.3.5 | sklearn 1.9.0 | xgboost 1.6.2 | lightgbm 4.7.0
XGBoost 3-fold CV ROC-AUC = 0.9042 +/- 0.0007  (mean +/- std across 3 stratified folds; a small std means a stable estimate on this split)

Reported approach	Known caveat
Wiefling et al. (2022) RBA scoring	IP/ASN reputation carries much of it
Freeman et al. (2016) authenticity model	features overlap with the labelling signal

Property	Value
Source	Kaggle `dasgroup/rba-dataset` (Wiefling et al., 2022)
Rows	31.3M logins; every attack row kept + benign subsampled to ≥1M
Label	`Is Attack IP` (the only target used); `Is Account Takeover` (loaded, then dropped, never modelled)
Access	Kaggle API token required (~9 GB one-time download)

	setting	held_out_auc
0	headline (as-is)	0.909410
1	de-duplicated (90% rows removed)	0.947699
2	shortcut feature dropped (ASN)	0.888466

Goal of this notebook¶

What you will learn¶

Where this connects to the course text¶

1. Research problem¶

2. Literature review¶

3. Dataset provenance & honesty caveats¶

Before you run this: getting the data¶

4. Solution design¶

5. Implementation architecture¶

6. Data acquisition & preparation¶

7. Exploratory data analysis¶

8. Model comparison¶

9. Results¶

10. Validity audit — is the score real?¶

11. Ablation — does the headline survive removing the artifacts?¶

12. Reproducibility & robustness¶

13. Scientific conclusion¶

References¶

	model	accuracy	roc_auc	train_s
0	XGBoost	0.859547	0.909410	0.5
1	LightGBM	0.857329	0.908066	1.5
2	RandomForest	0.859579	0.907558	1.3
3	LogisticRegression	0.786163	0.815961	0.3
4	MajorityBaseline	0.521900	0.500000	0.0

Goal of this notebook¶

What you will learn¶

Where this connects to the course text¶

Credential-Attack Detection on the RBA Login Dataset¶

Model comparison + validity audit on ≥1M login attempts (RBA synthesized release, 31.3M total)¶

1. Research problem¶

2. Literature review¶

3. Dataset provenance & honesty caveats¶

Before you run this: getting the data¶

4. Solution design¶

5. Implementation architecture¶

6. Data acquisition & preparation¶

7. Exploratory data analysis¶

8. Model comparison¶

9. Results¶

10. Validity audit — is the score real?¶

11. Ablation — does the headline survive removing the artifacts?¶

12. Reproducibility & robustness¶

13. Scientific conclusion¶

References¶