Data Science

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

df = pd.read_csv("data.csv")

# Shape, types, missing values
print(df.shape)
print(df.dtypes)
print(df.isnull().sum().sort_values(ascending=False))

# Numeric summary
print(df.describe())

# Categorical value counts
for col in df.select_dtypes("object"):
    print(f"\n{col}:\n{df[col].value_counts().head(10)}")

# Distribution of each numeric feature
df.hist(bins=30, figsize=(14, 10))
plt.tight_layout()
plt.show()

# Correlation heatmap
plt.figure(figsize=(10, 8))
sns.heatmap(
    df.select_dtypes("number").corr(),
    annot=True, fmt=".2f", cmap="coolwarm", center=0
)
plt.show()

def clean_dataframe(df: pd.DataFrame) -> pd.DataFrame:
    df = df.copy()  # Never mutate the original

    # 1. Standardize column names
    df.columns = df.columns.str.lower().str.replace(r"\s+", "_", regex=True)

    # 2. Drop duplicates
    df = df.drop_duplicates()

    # 3. Handle missing values
    numeric_cols = df.select_dtypes("number").columns
    categorical_cols = df.select_dtypes("object").columns

    df[numeric_cols] = df[numeric_cols].fillna(df[numeric_cols].median())
    df[categorical_cols] = df[categorical_cols].fillna("unknown")

    # 4. Remove outliers (IQR method - only for numeric targets)
    for col in numeric_cols:
        q1, q3 = df[col].quantile([0.25, 0.75])
        iqr = q3 - q1
        df = df[(df[col] >= q1 - 1.5 * iqr) & (df[col] <= q3 + 1.5 * iqr)]

    return df

from scipy import stats

# Independent samples t-test (two groups, continuous outcome)
group_a = df[df["variant"] == "control"]["revenue"]
group_b = df[df["variant"] == "treatment"]["revenue"]

# Check normality (Shapiro-Wilk - only reliable for n < 5000)
_, p_norm_a = stats.shapiro(group_a.sample(min(len(group_a), 500)))
_, p_norm_b = stats.shapiro(group_b.sample(min(len(group_b), 500)))
print(f"Normality p-values: A={p_norm_a:.4f}, B={p_norm_b:.4f}")

# If p_norm < 0.05 on small samples, prefer Mann-Whitney U
if p_norm_a < 0.05 or p_norm_b < 0.05:
    stat, p_value = stats.mannwhitneyu(group_a, group_b, alternative="two-sided")
    print(f"Mann-Whitney U: stat={stat:.2f}, p={p_value:.4f}")
else:
    stat, p_value = stats.ttest_ind(group_a, group_b)
    print(f"t-test: t={stat:.2f}, p={p_value:.4f}")

# Effect size (Cohen's d)
pooled_std = np.sqrt((group_a.std() ** 2 + group_b.std() ** 2) / 2)
cohens_d = (group_b.mean() - group_a.mean()) / pooled_std
print(f"Cohen's d: {cohens_d:.3f}")  # < 0.2 small, 0.5 medium, > 0.8 large

# Chi-square test for categorical outcomes
contingency = pd.crosstab(df["variant"], df["converted"])
chi2, p_chi2, dof, expected = stats.chi2_contingency(contingency)
print(f"Chi-square: chi2={chi2:.2f}, p={p_chi2:.4f}, dof={dof}")

from statsmodels.stats.power import TTestIndPower, NormalIndPower
from statsmodels.stats.proportion import proportions_ztest

# Sample size for conversion rate test
# effect_size = (p2 - p1) / sqrt(p_pooled * (1 - p_pooled))
baseline_rate = 0.05       # current conversion
minimum_detectable = 0.01  # smallest change worth detecting
alpha = 0.05               # false positive rate
power = 0.80               # 1 - false negative rate

p1, p2 = baseline_rate, baseline_rate + minimum_detectable
p_pool = (p1 + p2) / 2
effect_size = (p2 - p1) / np.sqrt(p_pool * (1 - p_pool))

analysis = NormalIndPower()
n = analysis.solve_power(effect_size=effect_size, alpha=alpha, power=power)
print(f"Required n per group: {int(np.ceil(n))}")

# Analysis after experiment
control_conversions = 520
control_n = 10000
treatment_conversions = 570
treatment_n = 10000

counts = np.array([treatment_conversions, control_conversions])
nobs = np.array([treatment_n, control_n])
z_stat, p_value = proportions_ztest(counts, nobs)
lift = (treatment_conversions / treatment_n) / (control_conversions / control_n) - 1
print(f"Lift: {lift:.1%}, z={z_stat:.2f}, p={p_value:.4f}")

import matplotlib.pyplot as plt
import seaborn as sns

# Set a consistent style once at the top of the notebook
sns.set_theme(style="whitegrid", palette="muted", font_scale=1.1)

# Distribution comparison - violin > box when showing distribution shape
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
sns.violinplot(data=df, x="group", y="value", ax=axes[0])
axes[0].set_title("Distribution by Group")

# Scatter with regression line - always show the uncertainty band
sns.regplot(data=df, x="feature", y="target", scatter_kws={"alpha": 0.3}, ax=axes[1])
axes[1].set_title("Feature vs Target")
plt.tight_layout()

# Time series - always label axes and use ISO date format
fig, ax = plt.subplots(figsize=(12, 4))
ax.plot(df["date"], df["metric"], color="steelblue", linewidth=1.5)
ax.fill_between(df["date"], df["lower_ci"], df["upper_ci"], alpha=0.2)
ax.set_xlabel("Date")
ax.set_ylabel("Metric")
ax.set_title("Metric Over Time with 95% CI")
plt.xticks(rotation=45)
plt.tight_layout()

from sklearn.preprocessing import StandardScaler, LabelEncoder
from sklearn.model_selection import train_test_split

# 1. Split first - to prevent leakage
X = df.drop("target", axis=1)
y = df["target"]
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42, stratify=y
)

# 2. Numeric features - fit scaler on train, transform both
scaler = StandardScaler()
num_cols = X_train.select_dtypes("number").columns
X_train[num_cols] = scaler.fit_transform(X_train[num_cols])
X_test[num_cols] = scaler.transform(X_test[num_cols])  # transform only, no fit

# 3. Date features
df["hour"] = pd.to_datetime(df["timestamp"]).dt.hour
df["day_of_week"] = pd.to_datetime(df["timestamp"]).dt.dayofweek
df["is_weekend"] = df["day_of_week"].isin([5, 6]).astype(int)

# 4. Interaction features (only when domain knowledge suggests it)
df["price_per_sqft"] = df["price"] / df["sqft"].replace(0, np.nan)

# 5. Target encoding (use cross-val folds to prevent leakage)
from category_encoders import TargetEncoder
encoder = TargetEncoder(smoothing=10)
X_train["cat_encoded"] = encoder.fit_transform(X_train["category"], y_train)
X_test["cat_encoded"] = encoder.transform(X_test["category"])

from sklearn.linear_model import LinearRegression, LogisticRegression
from sklearn.metrics import (
    mean_squared_error, r2_score,
    classification_report, roc_auc_score
)
import statsmodels.api as sm

# Linear regression with statistical output (p-values, CIs)
X_with_const = sm.add_constant(X_train[["feature_1", "feature_2"]])
ols_model = sm.OLS(y_train, X_with_const).fit()
print(ols_model.summary())  # Shows coefficients, p-values, R-squared

# Sklearn for prediction pipeline
lr = LinearRegression()
lr.fit(X_train[num_cols], y_train)
y_pred = lr.predict(X_test[num_cols])
print(f"RMSE: {mean_squared_error(y_test, y_pred, squared=False):.4f}")
print(f"R2: {r2_score(y_test, y_pred):.4f}")

# Logistic regression
clf = LogisticRegression(max_iter=1000, random_state=42)
clf.fit(X_train[num_cols], y_train)
y_prob = clf.predict_proba(X_test[num_cols])[:, 1]
print(classification_report(y_test, clf.predict(X_test[num_cols])))
print(f"ROC-AUC: {roc_auc_score(y_test, y_prob):.4f}")