350 KiB
350 KiB
In [81]:
# Import necessary libraries
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from sklearn.metrics import mean_squared_error
from sklearn.utils import resample
from numpy.linalg import LinAlgErrorIn [82]:
# Set a seed for reproducibility of random operations like train-test split
seed = 22980254
# Load the dataset into a pandas DataFrame
df = pd.read_csv('housing_data.csv')In [83]:
# Basic information about the dataset
data_info = df.info()
# Initial exploration of specified variables
exploration = df[['medv', 'rm', 'rad']].describe()
# Checking for missing values
missing_values = df[['medv', 'rm', 'rad']].isnull().sum()<class 'pandas.core.frame.DataFrame'> RangeIndex: 506 entries, 0 to 505 Data columns (total 12 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 crim 506 non-null float64 1 zn 506 non-null float64 2 indus 506 non-null float64 3 nox 506 non-null float64 4 rm 506 non-null float64 5 age 506 non-null float64 6 dis 506 non-null float64 7 rad 506 non-null int64 8 tax 506 non-null int64 9 ptratio 506 non-null float64 10 lstat 506 non-null float64 11 medv 506 non-null float64 dtypes: float64(10), int64(2) memory usage: 47.6 KB
In [84]:
# Histograms
plt.figure(figsize=(15, 5))
plt.subplot(1, 3, 1)
sns.histplot(df['medv'], kde=True)
plt.title('Distribution of MEDV')
plt.subplot(1, 3, 2)
sns.histplot(df['rm'], kde=True)
plt.title('Distribution of RM')
plt.subplot(1, 3, 3)
sns.histplot(df['rad'], kde=True, bins=24)
plt.title('Distribution of RAD')
plt.tight_layout()
plt.show()
In [85]:
# Scatter plots
plt.figure(figsize=(10, 5))
plt.subplot(1, 2, 1)
sns.scatterplot(x='rm', y='medv', data=df)
plt.title('MEDV vs. RM')
plt.subplot(1, 2, 2)
sns.scatterplot(x='rad', y='medv', data=df)
plt.title('MEDV vs. RAD')
plt.tight_layout()
plt.show()
In [86]:
# Generate a correlation matrix heatmap to visualize the correlations between all variables
plt.figure(figsize=(14, 10))
correlation_matrix = df[['medv', 'rm', 'rad']].corr()
sns.heatmap(correlation_matrix, annot=True, cmap='coolwarm')
plt.title('Correlation Matrix Heatmap')
print(correlation_matrix)
medv rm rad medv 1.000000 0.695360 -0.381626 rm 0.695360 1.000000 -0.209847 rad -0.381626 -0.209847 1.000000
In [87]:
print(exploration)
print(missing_values)
medv rm rad count 506.000000 506.000000 506.000000 mean 22.532806 6.284634 9.549407 std 9.197104 0.702617 8.707259 min 5.000000 3.561000 1.000000 25% 17.025000 5.885500 4.000000 50% 21.200000 6.208500 5.000000 75% 25.000000 6.623500 24.000000 max 50.000000 8.780000 24.000000 medv 0 rm 0 rad 0 dtype: int64
In [88]:
# Define a custom function to split the data into training and testing sets
# This function shuffles the indices and splits the data accordingly, ensuring reproducibility by setting a random seed
import numpy as np
import pandas as pd
def custom_split(data, target_column, test_size=0.1, student_id=seed):
# Validate inputs
if len(data) == 0 or test_size <= 0 or test_size >= 1:
raise ValueError("Invalid data size or test size")
if target_column not in data.columns:
raise ValueError("Target column not found in data")
# Set random seed using student ID for reproducibility
np.random.seed(student_id)
# Shuffle the dataset more thoroughly
data_shuffled = data.sample(frac=1, random_state=student_id).reset_index(drop=True)
# Determine the size of the test set
test_set_size = int(len(data) * test_size)
# Ensure test_set_size is not larger than the dataset
test_set_size = min(test_set_size, len(data) - 1)
# Split the indices for the test and training sets
test_indices = np.arange(test_set_size)
train_indices = np.arange(test_set_size, len(data))
# Split the data into training and test sets
test_set = data_shuffled.iloc[test_indices.tolist()]
train_set = data_shuffled.iloc[train_indices.tolist()]
# Split features and target variable
X_train = train_set.drop(columns=target_column)
y_train = train_set[target_column]
X_test = test_set.drop(columns=target_column)
y_test = test_set[target_column]
return X_train, X_test, y_train, y_test
In [89]:
# Custom function to fit a linear regression model using the normal equation
def custom_lr(X, y):
X_b = np.c_[np.ones((X.shape[0], 1)), X] # Adding a column of ones for the intercept term
theta_best = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y) # Calculating best fit parameters
return theta_best
# Custom function to predict using the linear regression model
def custom_predict(X, theta):
X_b = np.c_[np.ones((X.shape[0], 1)), X]
return X_b.dot(theta)
# Custom function to calculate Mean Squared Error
def custom_mean_squared_error(y_true, y_pred):
mse = np.mean((y_true - y_pred) ** 2)
return mseIn [90]:
# Splitting the data
X_train, X_test, y_train, y_test = custom_split(df, 'medv')
In [91]:
# Fitting the model using the custom function
theta_best = custom_lr(X_train, y_train)In [92]:
# Coefficient Interpretation
coefficients = pd.DataFrame([theta_best[1:]], columns=X_train.columns, index=['Coefficient']).T
In [93]:
# Bootstrap Analysis
def bootstrap_analysis(data, n_bootstrap=1000):
bootstrap_coefs = []
for _ in range(n_bootstrap):
sample_data = data.sample(n=len(data), replace=True)
X_sample, y_sample = sample_data.drop(columns='medv'), sample_data['medv']
theta_sample = custom_lr(X_sample, y_sample)
bootstrap_coefs.append(theta_sample[1:]) # Exclude the intercept
return np.array(bootstrap_coefs)
bootstrap_coefs = bootstrap_analysis(df)In [94]:
# Calculating statistics for bootstrap coefficients
coef_std = np.std(bootstrap_coefs, axis=0)
confidence_intervals = np.percentile(bootstrap_coefs, [2.5, 97.5], axis=0)
In [95]:
# Preparing results for display
coef_analysis = pd.DataFrame({
'Coefficient Mean': np.mean(bootstrap_coefs, axis=0),
'Std Dev': coef_std,
'95% CI Lower': confidence_intervals[0],
'95% CI Upper': confidence_intervals[1]
}, index=X_train.columns)In [96]:
# Check for significance
coef_analysis['Significant'] = (coef_analysis['95% CI Lower'] > 0) | (coef_analysis['95% CI Upper'] < 0)
In [97]:
# Predicting and calculating MSE
y_pred_custom = custom_predict(X_test, theta_best)
mse_custom = custom_mean_squared_error(y_test, y_pred_custom)In [98]:
# Displaying the results
print("Custom Coefficients:", coefficients)
print("Custom Mean Squared Error on Test Set:", mse_custom)
Custom Coefficients: Coefficient crim -0.131907 zn 0.048926 indus 0.024283 nox -18.026122 rm 3.520948 age 0.003638 dis -1.535510 rad 0.334010 tax -0.013442 ptratio -1.004608 lstat -0.583388 Custom Mean Squared Error on Test Set: 21.57614518279428
In [99]:
for i, feature in enumerate(X_train.columns):
print(f"{feature}: {confidence_intervals[:, i]}")crim: [-0.17961279 -0.05838793] zn: [0.01795932 0.07418449] indus: [-0.06823981 0.13497858] nox: [-25.5537351 -11.27688918] rm: [2.27561259 5.28224894] age: [-0.02491452 0.03740492] dis: [-1.93131905 -1.0892111 ] rad: [0.18975246 0.42646345] tax: [-0.01970752 -0.00916768] ptratio: [-1.20534711 -0.75727759] lstat: [-0.75317524 -0.37464871]
In [100]:
# Perform bootstrap evaluation
def custom_bootstrap_evaluation(X_train, y_train, X_test, y_test, n_iterations=1000):
mse_values = []
for _ in range(n_iterations):
boot_x, boot_y = resample(X_train, y_train)
theta_boot = custom_lr(boot_x, boot_y)
y_pred_boot = custom_predict(X_test, theta_boot)
mse_values.append(custom_mean_squared_error(y_test, y_pred_boot))
return np.mean(mse_values), np.std(mse_values)
bootstrap_mse_mean, bootstrap_mse_std = custom_bootstrap_evaluation(X_train, y_train, X_test, y_test)
print(f"Custom Model - Bootstrap MSE Mean: {bootstrap_mse_mean}, Std: {bootstrap_mse_std}")Custom Model - Bootstrap MSE Mean: 22.266506448229325, Std: 1.9577850287579381
In [101]:
# Splitting the data using custom_split
X_train, X_test, y_train, y_test = custom_split(df, 'medv')
# Fitting the original model with all features
theta_best = custom_lr(X_train, y_train)
# Predicting on the test set with the original model
y_pred_original = custom_predict(X_test, theta_best)
mse_original = custom_mean_squared_error(y_test, y_pred_original)
# Refining the model to include only 'rm' and 'rad'
X_train_refined = X_train[['crim', 'zn', 'nox', 'rm', 'dis', 'rad', 'tax', 'ptratio', 'lstat']]
X_test_refined = X_test[['crim', 'zn', 'nox', 'rm', 'dis', 'rad', 'tax', 'ptratio', 'lstat']]
theta_refined = custom_lr(X_train_refined, y_train)
# Predicting on the test set with the refined model
y_pred_refined = custom_predict(X_test_refined, theta_refined)
mse_refined = custom_mean_squared_error(y_test, y_pred_refined)
# Coefficients of the refined model
refined_coefficients = pd.DataFrame([theta_refined[1:]], columns=X_train_refined.columns, index=['Coefficient']).T
print("Original Model MSE:", mse_original)
print("Refined Model MSE:", mse_refined)
# Plotting the results
plt.figure(figsize=(10, 5))
plt.bar(['Original Model', 'Refined Model'], [mse_original, mse_refined], color=['blue', 'red'])
plt.ylabel('Mean Squared Error')
plt.title('Original vs Refined Model MSE Comparison')
plt.show()
plt.figure(figsize=(10, 5))
plt.bar(refined_coefficients.index, refined_coefficients['Coefficient'], color='blue')
plt.ylabel('Coefficient Value')
plt.title('Coefficients of the Refined Model')
plt.show()Original Model MSE: 21.57614518279428 Refined Model MSE: 21.698343879488906
In [102]:
# Splitting the data
X_train, X_test, y_train, y_test = custom_split(df, 'medv')
# Fitting the original model
theta_best = custom_lr(X_train, y_train)
# Predicting on the test set with the original model
y_pred_original = custom_predict(X_test, theta_best)
# Refining the model to include only 'rm' and 'rad'
X_train_refined = X_train[['rm', 'rad']]
X_test_refined = X_test[['rm', 'rad']]
theta_refined = custom_lr(X_train_refined, y_train)
# Predicting on the test set with the refined model
y_pred_refined = custom_predict(X_test_refined, theta_refined)
# Preparing data for violin plot
plot_data = pd.DataFrame({
'Actual Values': y_test,
'Original Model Predictions': y_pred_original,
'Refined Model Predictions': y_pred_refined
})
# Melting the DataFrame for use with seaborn
plot_data_melted = plot_data.melt(var_name='Group', value_name='MEDV')
# Plotting violin plots
plt.figure(figsize=(12, 6))
sns.violinplot(x='Group', y='MEDV', data=plot_data_melted)
plt.title('Comparison of Prediction Distributions')
plt.show()
In [103]:
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from xgboost import XGBRegressor
from sklearn.metrics import mean_squared_error
import matplotlib.pyplot as plt
# Assuming df is your DataFrame and 'medv' is the target column
# Splitting the data using custom_split
X_train_custom, X_test_custom, y_train_custom, y_test_custom = custom_split(df, 'medv')
# Splitting the data using sklearn's train_test_split
X_train_sklearn, X_test_sklearn, y_train_sklearn, y_test_sklearn = train_test_split(df.drop(columns='medv'), df['medv'], test_size=0.1, random_state=42)
# Initialize lists to store MSEs
mse_custom_split = []
mse_sklearn_split = []
# Custom Linear Regression Model
theta_custom = custom_lr(X_train_custom, y_train_custom)
y_pred_custom = custom_predict(X_test_custom, theta_custom)
mse_custom_split.append(custom_mean_squared_error(y_test_custom, y_pred_custom))
theta_sklearn = custom_lr(X_train_sklearn, y_train_sklearn)
y_pred_sklearn = custom_predict(X_test_sklearn, theta_sklearn)
mse_sklearn_split.append(custom_mean_squared_error(y_test_sklearn, y_pred_sklearn))
# Sklearn Linear Regression Model
lr_model_custom = LinearRegression()
lr_model_custom.fit(X_train_custom, y_train_custom)
y_pred_lr_custom = lr_model_custom.predict(X_test_custom)
mse_custom_split.append(mean_squared_error(y_test_custom, y_pred_lr_custom))
lr_model_sklearn = LinearRegression()
lr_model_sklearn.fit(X_train_sklearn, y_train_sklearn)
y_pred_lr_sklearn = lr_model_sklearn.predict(X_test_sklearn)
mse_sklearn_split.append(mean_squared_error(y_test_sklearn, y_pred_lr_sklearn))
# XGBoost Model
xgb_model_custom = XGBRegressor()
xgb_model_custom.fit(X_train_custom, y_train_custom)
y_pred_xgb_custom = xgb_model_custom.predict(X_test_custom)
mse_custom_split.append(mean_squared_error(y_test_custom, y_pred_xgb_custom))
xgb_model_sklearn = XGBRegressor()
xgb_model_sklearn.fit(X_train_sklearn, y_train_sklearn)
y_pred_xgb_sklearn = xgb_model_sklearn.predict(X_test_sklearn)
mse_sklearn_split.append(mean_squared_error(y_test_sklearn, y_pred_xgb_sklearn))
# Plotting the results
plt.figure(figsize=(12, 6))
model_names = ['Custom Linear Regression', 'Sklearn Linear Regression', 'XGBoost']
bar_width = 0.35
index = np.arange(len(model_names))
plt.bar(index, mse_custom_split, bar_width, label='Custom Split')
plt.bar(index + bar_width, mse_sklearn_split, bar_width, label='Sklearn Split')
plt.xlabel('Model')
plt.ylabel('Mean Squared Error')
plt.title('MSE Comparison by Split Method')
plt.xticks(index + bar_width / 2, model_names)
plt.legend()
plt.tight_layout()
plt.show()
