<jupyter_start><jupyter_text>Verifying Central Limit Theorem in regression<jupyter_code>import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
import warnings
warnings.filterwarnings('ignore')<jupyter_output><empty_output><jupyter_text>Synthesize the datasetCreate `1000` random integers between `0`, `100` for `X` and create `y` such that$$y = \beta_{0} + \beta_{1}X + \epsilon$$where$$\beta_{0} = 30 \ and \ \beta_{1} = 1.8 \ and \ \epsilon \ = \ standard \ normal \ error$$<jupyter_code>rand_1kx = np.random.randint(0,100,1000)
x_mean = np.mean(rand_1kx)
x_sd = np.std(rand_1kx)
x_mean
pop_intercept = 30
pop_slope = 1.8
error_boost = 10
pop_error = np.random.standard_normal(size = rand_1kx.size) * error_boost
# I added an error booster since without it, the correlation was too high.

y = pop_intercept + pop_slope*rand_1kx + pop_error
y_mean = np.mean(y)
y_sd = np.std(y)
y_mean<jupyter_output><empty_output><jupyter_text>Make a scatter plot of `X` and `y` variables.<jupyter_code>sns.jointplot(rand_1kx, y)<jupyter_output><empty_output><jupyter_text>`X` and `y` follow `uniform` distribution, but the error $\epsilon$ is generated from `standard normal distribution` with a boosting factor. Let us plot its histogram to verify the distribution<jupyter_code>sns.distplot(pop_error)<jupyter_output><empty_output><jupyter_text>Predict using populationLet us predict the coefficients and intercept when using the whole dataset. We will compare this approach with CLT approach of breaking into multiple subsets and averaging the coefficients and intercepts  Using whole population<jupyter_code>from sklearn.linear_model import LinearRegression
X_train_full = rand_1kx.reshape(-1,1)
y_train_full = y.reshape(-1,1)
y_train_full.shape
lm.fit(X_train, y_train)

#print the linear model built
predicted_pop_slope = lm.coef_[0][0]
predicted_pop_intercept = lm.intercept_[0]

print("y = " + str(predicted_pop_slope) + "*X" + " + " + str(predicted_pop_intercept))<jupyter_output>y = 1.795560991921382*X + 30.718916711669976<jupyter_text>Prediction with 66% of data<jupyter_code>from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(rand_1kx, y, test_size=0.33)
print(X_train.size)

from sklearn.linear_model import LinearRegression
lm = LinearRegression()
X_train = X_train.reshape(-1,1)
X_test = X_test.reshape(-1,1)
y_train = y_train.reshape(-1,1)
y_test = y_test.reshape(-1,1)
y_train.shape
lm.fit(X_train, y_train)

#print the linear model built
predicted_subset_slope = lm.coef_[0][0]
predicted_subset_intercept = lm.intercept_[0]

print("y = " + str(predicted_subset_slope) + "*X" 
      + " + " + str(predicted_subset_intercept))<jupyter_output>y = 1.794887898705644*X + 29.857924099881075<jupyter_text>Perform predictions and plot the charts<jupyter_code>y_predicted = lm.predict(X_test)
residuals = y_test - y_predicted<jupyter_output><empty_output><jupyter_text>Fitted vs Actual scatter<jupyter_code>jax = sns.jointplot(y_test, y_predicted)
jax.set_axis_labels(xlabel='Y', ylabel='Predicted Y')
dax = sns.distplot(residuals)
dax.set_title('Distribution of residuals')
jax = sns.jointplot(y_predicted, residuals)
jax.set_axis_labels(xlabel='Predicted Y', ylabel='Residuals')
jax = sns.jointplot(y_test, residuals)
jax.set_axis_labels(xlabel='Y', ylabel='Residuals')<jupyter_output><empty_output><jupyter_text>Predict using multiple samples<jupyter_code>pop_df = pd.DataFrame(data={'x':rand_1kx, 'y':y})
pop_df.head()
pop_df.shape<jupyter_output><empty_output><jupyter_text>Select 50 samples of size 200 and perform regression<jupyter_code>sample_slopes = []
sample_intercepts = []

for i in range(0,50):
    # perform a choice on dataframe index
    sample_index = np.random.choice(pop_df.index, size=50)
    
    # select the subset using that index
    sample_df = pop_df.iloc[sample_index]
    
    # convert to numpy and reshape the matrix for lm.fit
    sample_x = np.array(sample_df['x']).reshape(-1,1)
    sample_y = np.array(sample_df['y']).reshape(-1,1)
    
    lm.fit(X=sample_x, y=sample_y)
    
    sample_slopes.append(lm.coef_[0][0])
    sample_intercepts.append(lm.intercept_[0])<jupyter_output><empty_output><jupyter_text>Plot the distribution of sample slopes and intercepts<jupyter_code>mean_sample_slope = np.mean(sample_slopes)
mean_sample_intercept = np.mean(sample_intercepts)

fig, ax = plt.subplots(1,2, figsize=(15,6))

# plot sample slopes
sns.distplot(sample_slopes, ax=ax[0])
ax[0].set_title('Distribution of sample slopes. Mean: ' 
                + str(round(mean_sample_slope, 2)))
ax[0].axvline(mean_sample_slope, color='black')

# plot sample slopes
sns.distplot(sample_intercepts, ax=ax[1])
ax[1].set_title('Distribution of sample intercepts. Mean: ' 
                + str(round(mean_sample_intercept,2)))
ax[1].axvline(mean_sample_intercept, color='black')<jupyter_output><empty_output><jupyter_text>ConclusionHere we compare the coefficients and intercepts obtained by different methods to see how CLT adds up.<jupyter_code>print("Predicting using population")
print("----------------------------")
print("Error in intercept: {}".format(pop_intercept - predicted_pop_intercept))
print("Error in slope: {}".format(pop_slope - predicted_pop_slope))

print("\n\nPredicting using subset")
print("----------------------------")
print("Error in intercept: {}".format(pop_intercept - predicted_subset_intercept))
print("Error in slope: {}".format(pop_slope - predicted_subset_slope))

print("\n\nPredicting using a number of smaller samples")
print("------------------------------------------------")
print("Error in intercept: {}".format(pop_intercept - mean_sample_intercept))
print("Error in slope: {}".format(pop_slope - mean_sample_slope))<jupyter_output>Predicting using population
----------------------------
Error in intercept: -0.7189167116699764
Error in slope: 0.0044390080786180786


Predicting using subset
----------------------------
Error in intercept: 0.14207590011892535
Error in slope: 0.0051121012943560196


Predicting using a number of smaller samples
------------------------------------------------
Error in intercept: 0.4823977050074646
Error in slope: 0.002971759530004725
{'fraction_non_alphanumeric': 0.14418455623197693, 'fraction_numerical': 0.043095161807113105, 'mean_word_length': 5.1628825271470875, 'pattern_counts': {'":': 0, '<': 46, '<?xml version=': 0, '>': 46, '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': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '133824', 'license': 'MIT', 'n_tokens_mistral': 2360, 'n_tokens_neox': 2183, 'n_words': 460, 'path': 'islr/verifying_clt_in_regression.ipynb', 'repo_name': 'AtmaMani/pyChakras', 'size': '6242'}