<jupyter_start><jupyter_text>Delete-1 Jackknife Explained  Load packages<jupyter_code>import numpy as np
from sklearn.model_selection import train_test_split

import matplotlib.pyplot as plt
%matplotlib inline

import scipy.stats<jupyter_output><empty_output><jupyter_text>Load dataset<jupyter_code>from sklearn.datasets import load_boston
tmp = load_boston()
y = tmp.target
X = tmp.data[:,[5,12]]<jupyter_output><empty_output><jupyter_text>Data SplittingLet's reserve 50% for testing lateron<jupyter_code>X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.5, random_state=42)

N = len(y_train)<jupyter_output><empty_output><jupyter_text>Model SpecificationWe will use a linear regression model but it could be any model with a fixed set of parameters/weights to be estimated.<jupyter_code>def estimator_fn(y, X, bias=True):
    if bias:
        X_ = np.c_[np.ones(shape=(len(y), 1)), X]
    else:
        X_ = X
    # OLS by LU decomposition
    return np.linalg.solve(np.dot(X_.T, X_), np.dot(X_.T, y))
def predictor_fn(X, weights, bias=True):
    if bias:
        X_ = np.c_[np.ones(shape=(len(X), 1)), X]
    else:
        X_ = X
    return np.dot(X_, weights)
def loss_fn(y, y_pred):
    return np.mean((y - y_pred)**2)<jupyter_output><empty_output><jupyter_text>Model Estimation Loop over the SubsamplesIn Delete-1 Jackknife `N` subsamples are created that leaves out 1 of the `N` observations. Then the parameters/weights of the base model are estimated for each of the `N` subsamples.<jupyter_code>%%time
params = []  # store estimated weights lateron

# row indicies for each observation
idx = np.arange(0, N)

#loop over subsamples
for i in idx:
    # create true/false mask for the i-th subsample
    # there is one observation masked false!
    mask = idx != i

    #estimate the the model parameters for the i-th subsample
    beta = estimator_fn(y_train[mask], X_train[mask, :])
    params.append(beta)  # store it

params = np.vstack(params)  # convert to numpy array<jupyter_output>CPU times: user 42.3 ms, sys: 2.76 ms, total: 45.1 ms
Wall time: 89.1 ms<jupyter_text>Eyeball the subsamples' estimated parameters/weights ('Jackknife replicates')<jupyter_code>fig = plt.figure(figsize=(5,3), dpi=150)
names = ['intercept', 'beta_1', 'beta_2']
color = ['b+', 'r>', 'g*']
bwidth = .7

for i in range(3):
    ax = plt.subplot('13'+str(i+1))
    ax.boxplot(x=params[:,i], sym=color[i], widths=bwidth);
    ax.set_xlabel(names[i])
    ax.set_xticks([])

plt.tight_layout()<jupyter_output><empty_output><jupyter_text>Remember, the classic Delete-1 Jackknife just leaves out 1 observation for each subsample. This is already enough to create subsamples with fluctating parameter estimates (The Jackknife Replicates). For example, the *intercept* parameter estimates are centered around `-4.5` but one subsample resulted in an almost `-8.0` parameter value, or almost twice as negative. The Jackknife method literally *average out* the outliers.   Aggregate EstimatesIn a first step, the subsamples parameters/weights are aggregated:- $\theta_{biased}$ : Biased Jackknife Estimate (mean of weights)- $SE(\theta)$ : Jackknife Standard Error (standard deviation of weights)- $SE(\theta) / |\theta_{biased}|$ : Quick-n-Dirty indicator for parameter/weight stabilityAnd finally, the mathematically correct answer- $\theta_{jack}$ : Bias-corrected Jackknife Estimate (use this as model parameters/weights!)- Jackknife t-Statistics (reject the whole model if any parameter/weight is unstable)  Biased Jackknife EstimateThe biased **Jackknife Estimate** is the mean of a specific parameter estimate (e.g. the intercept) for all subsamples.$$\theta_{biased} = \frac{1}{n} \sum_{i=1}^n \theta_{(i)}$$In our example, the biased Jackknife Estimate is<jupyter_code>theta_biased = np.mean(params, axis=0)
theta_biased<jupyter_output><empty_output><jupyter_text>Jackknife Standard Error The core of the Jackknife procedure is the question how stable model parameters actually are. The **Standard Error (SE) of the Jackknife Estimates** is$$SE(\theta) = \sqrt{ \frac{n - 1}{n} \sum_{i=1}^n \left(\theta_{(i)} - \theta_{biased}\right)^2 }$$<jupyter_code>se_jack = np.sqrt( (N-1)/N * np.sum((params - theta_biased)**2, axis=0) )
se_jack<jupyter_output><empty_output><jupyter_text>Quick-n-Dirty indicator for parameter/weight stabilityWhat parameters are the most unstable?It is not really scientific but we can already guess potentially unstable parameter estimates without bias-correction by$$\frac{SE(\theta)}{|\theta_{biased}|}$$A low standard error per absolute parameter value is desirable (lower is better or more stable). In our example, the intercept seems problematic.<jupyter_code>se_jack / abs(theta_biased)<jupyter_output><empty_output><jupyter_text>Bias-corrected Jackknife EstimateThe **unbiased** or **bias-corrected Jackknife Estimate** is the parameter estimate of the full sample minus a Jackknife specific bias-correction$$\theta_{jack} = \theta_{all} - bias_{jack}$$First, let estimate the model on the full sample<jupyter_code>theta_all = estimator_fn(y_train, X_train)
theta_all<jupyter_output><empty_output><jupyter_text>The **Jackknife Bias-Correction** $bias_{jack}$ is$$bias_{jack} = (n - 1) (\theta_{all} - \theta_{biased})$$<jupyter_code>bias_jack = (N - 1) * (theta_biased - theta_all)
bias_jack<jupyter_output><empty_output><jupyter_text>thus<jupyter_code>theta_jack = theta_all - bias_jack
theta_jack<jupyter_output><empty_output><jupyter_text>Substitution merge both formulas into$$\theta_{jack} = n \, \theta_{all} - (n - 1) \, \theta_{biased}$$In our example,<jupyter_code>theta_jack = N * theta_all - (N - 1) * theta_biased
theta_jack<jupyter_output><empty_output><jupyter_text>Jackknife t-Statistics - What parameters are unstable? The t-Statistics is simply<jupyter_code>tscore = theta_jack / se_jack
tscore<jupyter_output><empty_output><jupyter_text>The two-sided p-value is<jupyter_code>M = len(tscore);
pval = scipy.stats.t.sf(np.abs(tscore), N-M-1)*2
pval.round(4)<jupyter_output><empty_output><jupyter_text>And as we already guessed, the intercept parameter estimates are just random (`0.524` is way of off being close to Zero), i.e. the intercept parameter estimates in the Jackknife subsamples just happenend by *chance*.**As a consequence the whole model is rejected.**The analysis is done.   Prediction However, we could still burn through our test set. It does not matter anyways ;)<jupyter_code># theta_all estimated from full training set
y_pred  = predictor_fn(X_train, theta_all)
mse_all_in = loss_fn(y_train, y_pred)

y_pred  = predictor_fn(X_test, theta_all)
mse_all_out = loss_fn(y_test, y_pred)

mse_all_in, mse_all_out
# theta_jack, the biased corrected jackknife weights
y_pred  = predictor_fn(X_train, theta_jack)
mse_jack_in = loss_fn(y_train, y_pred)

y_pred  = predictor_fn(X_test, theta_jack)
mse_jack_out = loss_fn(y_test, y_pred)

mse_jack_in, mse_jack_out<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11272257323377369, 'fraction_numerical': 0.009477311889718553, 'mean_word_length': 5.749031007751938, 'pattern_counts': {'":': 0, '<': 64, '<?xml version=': 0, '>': 65, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '344211', 'license': 'MIT', 'n_tokens_mistral': 2432, 'n_tokens_neox': 2317, 'n_words': 667, 'path': 'examples/Delete-1 Jackknife Explained.ipynb', 'repo_name': 'kmedian/jackknife', 'size': '6964'}