Classification with Scikit-Learn — Algorithms, Pipelines, and Evaluation

from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import classification_report

X, y = load_iris(return_X_y=True)

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42, stratify=y
)

clf = RandomForestClassifier(n_estimators=100, random_state=42)
clf.fit(X_train, y_train)

y_pred = clf.predict(X_test)
print(classification_report(y_test, y_pred, target_names=['setosa', 'versicolor', 'virginica']))

from sklearn.linear_model import LogisticRegression
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import RandomForestClassifier, GradientBoostingClassifier
from sklearn.svm import SVC

classifiers = {
    'Logistic Regression': LogisticRegression(max_iter=1000),
    'Decision Tree':       DecisionTreeClassifier(max_depth=5),
    'Random Forest':       RandomForestClassifier(n_estimators=100),
    'Gradient Boosting':   GradientBoostingClassifier(n_estimators=100),
    'SVM':                 SVC(probability=True),
}

for name, clf in classifiers.items():
    clf.fit(X_train, y_train)
    acc = clf.score(X_test, y_test)
    print(f'{name}: {acc:.4f}')

probs = RandomForestClassifier().fit(X_train, y_train).predict_proba(X_test)
print(probs[0])

from sklearn.preprocessing import StandardScaler, MinMaxScaler, RobustScaler
from sklearn.pipeline import Pipeline
from sklearn.svm import SVC
from sklearn.neural_network import MLPClassifier

# SVM NEEDS scaling
svm_scaled = Pipeline([
    ('scaler', StandardScaler()),
    ('clf', SVC(kernel='rbf', probability=True)),
])

# When data has outliers, use RobustScaler
svm_robust = Pipeline([
    ('scaler', RobustScaler()),
    ('clf', SVC(kernel='rbf', probability=True)),
])

# For neural networks or bounded-input models
nn_scaled = Pipeline([
    ('scaler', MinMaxScaler()),
    ('clf', MLPClassifier(hidden_layer_sizes=(100, 50), max_iter=500)),
])

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler, OneHotEncoder, OrdinalEncoder
from sklearn.impute import SimpleImputer
from sklearn.compose import ColumnTransformer
from sklearn.ensemble import RandomForestClassifier
import numpy as np
import pandas as pd

# Simulate a realistic dataset
df = pd.DataFrame({
    'age': [25, 45, np.nan, 30, 60],
    'income': [50000, 80000, 120000, np.nan, 90000],
    'country': ['US', 'UK', 'US', 'DE', 'FR'],
    'plan_type': ['basic', 'premium', 'basic', 'enterprise', 'premium'],
    'education': ['high_school', 'bachelors', 'masters', 'phd', 'bachelors'],
    'churned': [0, 0, 1, 0, 1],
})

X = df.drop('churned', axis=1)
y = df['churned']

numeric_features = ['age', 'income']
nominal_features = ['country', 'plan_type']
ordinal_features = ['education']

# Ordinal encoding for ordered categories
education_order = ['high_school', 'bachelors', 'masters', 'phd']

preprocessor = ColumnTransformer([
    ('num', Pipeline([('imputer', SimpleImputer(strategy='median')),
                      ('scaler', StandardScaler())]), numeric_features),
    ('nom', Pipeline([('imputer', SimpleImputer(strategy='most_frequent')),
                      ('onehot', OneHotEncoder(handle_unknown='ignore'))]), nominal_features),
    ('ord', Pipeline([('imputer', SimpleImputer(strategy='most_frequent')),
                      ('ordinal', OrdinalEncoder(categories=[education_order]))]), ordinal_features),
])

pipeline = Pipeline([
    ('preprocessor', preprocessor),
    ('classifier', RandomForestClassifier(n_estimators=100, random_state=42)),
])

pipeline.fit(X, y)

# See transformed feature names
feature_names = pipeline.named_steps['preprocessor'].get_feature_names_out()
print('Transformed features:', feature_names)

from sklearn.naive_bayes import GaussianNB, MultinomialNB, BernoulliNB
from sklearn.feature_extraction.text import TfidfVectorizer
from sklearn.pipeline import Pipeline
from sklearn.metrics import classification_report

# Gaussian Naive Bayes for continuous features
gnb = GaussianNB()
gnb.fit(X_train, y_train)
y_pred = gnb.predict(X_test)
print('GaussianNB:', classification_report(y_test, y_pred))

# Multinomial Naive Bayes for text (TF-IDF features)
text_pipeline = Pipeline([
    ('tfidf', TfidfVectorizer(max_features=5000, stop_words='english')),
    ('clf',   MultinomialNB(alpha=0.1)),
])

# Bernoulli Naive Bayes for binary features
bnb = BernoulliNB()
binary_features = (X_train > 0).astype(int)
bnb.fit(binary_features, y_train)

from sklearn.metrics import (
    classification_report, confusion_matrix,
    roc_auc_score, average_precision_score,
    ConfusionMatrixDisplay, precision_recall_curve
)
import matplotlib.pyplot as plt
import numpy as np

print(classification_report(y_test, y_pred))

cm = confusion_matrix(y_test, y_pred)
disp = ConfusionMatrixDisplay(confusion_matrix=cm)
disp.plot()
plt.title('Confusion Matrix')
plt.show()

# ROC-AUC for binary classification
probs = pipeline.predict_proba(X_test)[:, 1]
auc = roc_auc_score(y_test, probs)
print(f'ROC-AUC: {auc:.4f}')

# PR-AUC — better metric for imbalanced classes
pr_auc = average_precision_score(y_test, probs)
print(f'PR-AUC:  {pr_auc:.4f}')

# Custom class weights — when you know the business cost ratio
clf_cost_sensitive = RandomForestClassifier(
    class_weight={0: 1, 1: 100},
    random_state=42
)

from sklearn.metrics import precision_recall_curve, f1_score
import numpy as np

y_probs = pipeline.predict_proba(X_test)[:, 1]
precisions, recalls, thresholds = precision_recall_curve(y_test, y_probs)

# Find threshold that maximises F1
f1_scores = 2 * (precisions * recalls) / (precisions + recalls + 1e-8)
best_idx = np.argmax(f1_scores)
best_threshold = thresholds[best_idx]
print(f'Best threshold for F1: {best_threshold:.3f}')
print(f'Precision: {precisions[best_idx]:.3f}, Recall: {recalls[best_idx]:.3f}, F1: {f1_scores[best_idx]:.3f}')

# Apply custom threshold
y_pred_custom = (y_probs >= best_threshold).astype(int)
print(f'F1 with custom threshold: {f1_score(y_test, y_pred_custom):.4f}')

# Find threshold for minimum recall target
min_recall = 0.95
valid_indices = np.where(recalls[:-1] >= min_recall)[0]
if len(valid_indices) > 0:
    best_for_recall = valid_indices[np.argmax(precisions[valid_indices])]
    recall_threshold = thresholds[best_for_recall]
    print(f'Threshold for >= 95% recall: {recall_threshold:.3f}')

# Business cost minimisation
fn_cost = 1000  # missed fraud
fp_cost = 10    # false alarm
min_cost = float('inf')
best_cost_threshold = 0.5
for i, t in enumerate(thresholds):
    y_pred_t = (y_probs >= t).astype(int)
    fn = np.sum((y_test == 1) & (y_pred_t == 0))
    fp = np.sum((y_test == 0) & (y_pred_t == 1))
    cost = fn * fn_cost + fp * fp_cost
    if cost < min_cost:
        min_cost = cost
        best_cost_threshold = t
print(f'Threshold minimising business cost: {best_cost_threshold:.3f} (cost: ${min_cost:,.0f})')

from imblearn.over_sampling import SMOTE, ADASYN
from imblearn.under_sampling import RandomUnderSampler
from imblearn.combine import SMOTEENN
from imblearn.pipeline import Pipeline as ImbPipeline
from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import StandardScaler, OneHotEncoder
from sklearn.impute import SimpleImputer
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import cross_val_score, StratifiedKFold
import numpy as np

# Strategy 1: class_weight (simplest)
clf_weighted = RandomForestClassifier(n_estimators=200, class_weight='balanced', random_state=42)

# Strategy 2: SMOTE inside imblearn pipeline
preprocessor = ColumnTransformer([
    ('num', Pipeline([('imp', SimpleImputer(strategy='median')), ('sc', StandardScaler())]), numeric_features),
    ('cat', Pipeline([('imp', SimpleImputer(strategy='most_frequent')), ('ohe', OneHotEncoder(handle_unknown='ignore'))]), categorical_features),
])

smote_pipeline = ImbPipeline([
    ('preprocessor', preprocessor),
    ('smote', SMOTE(random_state=42)),
    ('classifier', RandomForestClassifier(n_estimators=200, random_state=42)),
])

# Strategy 3: SMOTE + Edited Nearest Neighbours (cleans noisy samples)
smoteenn_pipeline = ImbPipeline([
    ('preprocessor', preprocessor),
    ('smoteenn', SMOTEENN(random_state=42)),
    ('classifier', RandomForestClassifier(n_estimators=200, random_state=42)),
])

# Strategy 4: Undersampling + SMOTE
combined_pipeline = ImbPipeline([
    ('preprocessor', preprocessor),
    ('under', RandomUnderSampler(sampling_strategy=0.5, random_state=42)),
    ('smote', SMOTE(random_state=42)),
    ('classifier', RandomForestClassifier(n_estimators=200, random_state=42)),
])

cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)

for name, pipe in [('class_weight', clf_weighted),
                   ('SMOTE', smote_pipeline),
                   ('SMOTEENN', smoteenn_pipeline),
                   ('Under+SMOTE', combined_pipeline)]:
    scores = cross_val_score(pipe, X_train, y_train, cv=cv, scoring='f1')
    print(f'{name:15s} F1: {scores.mean():.4f} (+/- {scores.std():.4f})')

from sklearn.inspection import permutation_importance
from sklearn.ensemble import RandomForestClassifier
from sklearn.feature_selection import SelectFromModel, SelectKBest, f_classif
from sklearn.pipeline import Pipeline
import pandas as pd
import numpy as np

clf = RandomForestClassifier(n_estimators=200, random_state=42)
clf.fit(X_train, y_train)

# Built-in feature importance
if hasattr(X_train, 'columns'):
    feature_names = X_train.columns
else:
    feature_names = [f'feature_{i}' for i in range(X_train.shape[1])]

importance_df = pd.DataFrame({
    'feature': feature_names,
    'importance': clf.feature_importances_
}).sort_values('importance', ascending=False)
print(importance_df.head(10))

# Permutation importance — more reliable
perm_imp = permutation_importance(clf, X_test, y_test, n_repeats=10, random_state=42, scoring='f1_weighted')
perm_df = pd.DataFrame({
    'feature': feature_names,
    'importance_mean': perm_imp.importances_mean,
    'importance_std': perm_imp.importances_std,
}).sort_values('importance_mean', ascending=False)
print(perm_df.head(10))

# Feature selection with SelectFromModel
selector = SelectFromModel(clf, threshold='median')
selector.fit(X_train, y_train)
X_selected = selector.transform(X_train)
print(f'Original features: {X_train.shape[1]}, Selected: {X_selected.shape[1]}')

# Pipeline with feature selection
pipeline_with_selection = Pipeline([
    ('selector', SelectFromModel(RandomForestClassifier(50, random_state=42), threshold='median')),
    ('classifier', RandomForestClassifier(200, random_state=42)),
])

import shap
from sklearn.ensemble import RandomForestClassifier

clf = RandomForestClassifier(n_estimators=200, random_state=42)
clf.fit(X_train, y_train)

# TreeExplainer — fast for tree-based models
explainer = shap.TreeExplainer(clf)
shap_values = explainer.shap_values(X_test)

# Global feature importance (mean absolute SHAP values)
shap.summary_plot(shap_values, X_test, feature_names=feature_names)

# Local explanation for a single prediction
sample_idx = 0
shap.force_plot(explainer.expected_value, shap_values[sample_idx],
                X_test.iloc[sample_idx], feature_names=feature_names)

# For multi-class, shap_values is a list of arrays
if isinstance(shap_values, list):
    print(f'SHAP values for class 0: shape {shap_values[0].shape}')
    print(f'SHAP values for class 1: shape {shap_values[1].shape}')

# Dependence plot: how does one feature affect predictions?
shap.dependence_plot('age', shap_values, X_test, feature_names=feature_names)

from sklearn.model_selection import cross_val_score, cross_validate, StratifiedKFold
from sklearn.ensemble import RandomForestClassifier

cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
clf = RandomForestClassifier(n_estimators=100, random_state=42)

scores = cross_val_score(clf, X, y, cv=cv, scoring='f1_weighted')
print(f'F1: {scores.mean():.4f} (+/- {scores.std():.4f})')

results = cross_validate(clf, X, y, cv=cv, scoring=['accuracy', 'f1_weighted'], return_train_score=True)
print('Train F1:', results['train_f1_weighted'].mean())
print('Test  F1:', results['test_f1_weighted'].mean())

# CV with full pipeline — no leakage
cv_scores = cross_val_score(pipeline, X, y, cv=cv, scoring='roc_auc')
print(f'Pipeline AUC: {cv_scores.mean():.4f}')

from sklearn.model_selection import GridSearchCV, RandomizedSearchCV, StratifiedKFold
from sklearn.ensemble import RandomForestClassifier
from scipy.stats import randint

cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)

# Stage 1: Random search
param_distributions = {
    'classifier__n_estimators': randint(50, 500),
    'classifier__max_depth': [None, 5, 10, 20, 30],
    'classifier__max_features': ['sqrt', 'log2', None],
    'classifier__min_samples_split': randint(2, 20),
}

random_search = RandomizedSearchCV(
    pipeline, param_distributions, n_iter=50,
    cv=cv, scoring='f1_weighted', n_jobs=-1, random_state=42
)
random_search.fit(X_train, y_train)
print('Best random params:', random_search.best_params_)
print('Best CV F1:', random_search.best_score_)

# Stage 2: Fine-grained grid around best values
best = random_search.best_params_
param_grid = {
    'classifier__n_estimators': [best['classifier__n_estimators'] - 50, best['classifier__n_estimators'], best['classifier__n_estimators'] + 50],
    'classifier__max_depth': [best['classifier__max_depth']],
}

grid_search = GridSearchCV(pipeline, param_grid, cv=cv, scoring='f1_weighted', n_jobs=-1)
grid_search.fit(X_train, y_train)
print('Best grid params:', grid_search.best_params_)

# Optuna
import optuna
from sklearn.model_selection import cross_val_score

def objective(trial):
    params = {
        'n_estimators': trial.suggest_int('n_estimators', 50, 500),
        'max_depth': trial.suggest_int('max_depth', 3, 30),
        'max_features': trial.suggest_categorical('max_features', ['sqrt', 'log2']),
        'min_samples_split': trial.suggest_int('min_samples_split', 2, 20),
    }
    clf = RandomForestClassifier(**params, random_state=42)
    return cross_val_score(clf, X_train, y_train, cv=5, scoring='f1_weighted').mean()

study = optuna.create_study(direction='maximize')
study.optimize(objective, n_trials=100)
print('Best Optuna params:', study.best_params)

from xgboost import XGBClassifier
from lightgbm import LGBMClassifier
from catboost import CatBoostClassifier
from sklearn.model_selection import cross_val_score, StratifiedKFold

cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)

xgb_clf = XGBClassifier(n_estimators=200, max_depth=6, learning_rate=0.1,
    subsample=0.8, colsample_bytree=0.8, reg_alpha=0.1, reg_lambda=1.0,
    random_state=42, eval_metric='logloss', use_label_encoder=False)

lgbm_clf = LGBMClassifier(n_estimators=200, max_depth=-1, learning_rate=0.1,
    num_leaves=31, subsample=0.8, reg_alpha=0.1, reg_lambda=1.0,
    random_state=42, verbose=-1)

cat_clf = CatBoostClassifier(iterations=200, depth=6, learning_rate=0.1,
    l2_leaf_reg=3.0, random_state=42, verbose=0)

for name, clf in [('XGBoost', xgb_clf), ('LightGBM', lgbm_clf), ('CatBoost', cat_clf)]:
    scores = cross_val_score(clf, X_train, y_train, cv=cv, scoring='f1_weighted')
    print(f'{name:10s} F1: {scores.mean():.4f} (+/- {scores.std():.4f})')

# CatBoost with native categorical features — no OneHotEncoder needed
cat_native = CatBoostClassifier(iterations=100, cat_features=['country', 'plan'], verbose=0)
cat_native.fit(X_train, y_train)

from sklearn.ensemble import VotingClassifier, StackingClassifier, BaggingClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.svm import SVC
from sklearn.tree import DecisionTreeClassifier
from sklearn.model_selection import cross_val_score, StratifiedKFold

cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)

voting_clf = VotingClassifier(
    estimators=[
        ('rf', RandomForestClassifier(100, random_state=42)),
        ('gb', GradientBoostingClassifier(100, random_state=42)),
        ('lr', LogisticRegression(max_iter=1000)),
    ],
    voting='soft',
)

stacking_clf = StackingClassifier(
    estimators=[
        ('rf', RandomForestClassifier(100, random_state=42)),
        ('gb', GradientBoostingClassifier(100, random_state=42)),
        ('svc', SVC(probability=True, random_state=42)),
    ],
    final_estimator=LogisticRegression(max_iter=1000),
    cv=5,
)

bagging_clf = BaggingClassifier(
    estimator=DecisionTreeClassifier(),
    n_estimators=100, random_state=42,
)

for name, clf in [('Voting', voting_clf), ('Stacking', stacking_clf), ('Bagging', bagging_clf)]:
    scores = cross_val_score(clf, X_train, y_train, cv=cv, scoring='f1_weighted')
    print(f'{name:10s} F1: {scores.mean():.4f} (+/- {scores.std():.4f})')

from sklearn.model_selection import learning_curve, validation_curve
from sklearn.ensemble import RandomForestClassifier
import matplotlib.pyplot as plt
import numpy as np

clf = RandomForestClassifier(n_estimators=100, random_state=42)
train_sizes, train_scores, val_scores = learning_curve(
    clf, X_train, y_train, train_sizes=np.linspace(0.1, 1.0, 10),
    cv=5, scoring='f1_weighted', n_jobs=-1)

fig, axes = plt.subplots(1, 2, figsize=(12, 5))
axes[0].plot(train_sizes, train_scores.mean(axis=1), label='Train')
axes[0].fill_between(train_sizes,
    train_scores.mean(axis=1) - train_scores.std(axis=1),
    train_scores.mean(axis=1) + train_scores.std(axis=1), alpha=0.1)
axes[0].plot(train_sizes, val_scores.mean(axis=1), label='Validation')
axes[0].fill_between(train_sizes,
    val_scores.mean(axis=1) - val_scores.std(axis=1),
    val_scores.mean(axis=1) + val_scores.std(axis=1), alpha=0.1)
axes[0].set_xlabel('Training Set Size')
axes[0].set_ylabel('F1 Score')
axes[0].set_title('Learning Curve')
axes[0].legend()

param_range = [2, 5, 10, 15, 20, 30]
train_vc, val_vc = validation_curve(
    RandomForestClassifier(100, random_state=42), X_train, y_train,
    param_name='max_depth', param_range=param_range, cv=5,
    scoring='f1_weighted', n_jobs=-1)

axes[1].plot(param_range, train_vc.mean(axis=1), label='Train')
axes[1].plot(param_range, val_vc.mean(axis=1), label='Validation')
axes[1].set_xlabel('max_depth')
axes[1].set_ylabel('F1 Score')
axes[1].set_title('Validation Curve')
axes[1].legend()
plt.tight_layout()
plt.show()

from sklearn.feature_extraction.text import TfidfVectorizer, CountVectorizer
from sklearn.naive_bayes import MultinomialNB
from sklearn.svm import LinearSVC
from sklearn.pipeline import Pipeline
from sklearn.model_selection import cross_val_score, StratifiedKFold

texts = [
    'Win a free iPhone now click here',
    'Meeting rescheduled to 3pm tomorrow',
    'Congratulations you won the lottery',
    'Please review the quarterly report',
    'Free pills online pharmacy discount',
    'Team standup notes for sprint 42',
    'Limited time offer act now',
    'Can you send the updated budget?',
    'You have been selected for a prize',
    'Project deadline extended to Friday',
] * 50
labels = [1, 0, 1, 0, 1, 0, 1, 0, 1, 0] * 50

# Naive Bayes + TF-IDF baseline
nb_pipeline = Pipeline([
    ('tfidf', TfidfVectorizer(max_features=5000, ngram_range=(1, 2), stop_words='english')),
    ('clf',   MultinomialNB(alpha=0.1)),
])

# SVM + TF-IDF — often better for text
svm_pipeline = Pipeline([
    ('tfidf', TfidfVectorizer(max_features=10000, ngram_range=(1, 2), stop_words='english')),
    ('clf',   LinearSVC(C=1.0, class_weight='balanced')),
])

cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)

for name, pipe in [('Naive Bayes', nb_pipeline), ('LinearSVC', svm_pipeline)]:
    scores = cross_val_score(pipe, texts, labels, cv=cv, scoring='f1')
    print(f'{name:12s} F1: {scores.mean():.4f} (+/- {scores.std():.4f})')

# Inspect most informative features
nb_pipeline.fit(texts, labels)
feature_names = nb_pipeline['tfidf'].get_feature_names_out()
for i, class_label in enumerate(['not spam', 'spam']):
    top_indices = nb_pipeline['clf'].feature_log_prob_[i].argsort()[-10:][::-1]
    print(f'Top words for {class_label}: {[feature_names[j] for j in top_indices]}')

from sklearn.calibration import CalibratedClassifierCV, calibration_curve
from sklearn.ensemble import RandomForestClassifier
import matplotlib.pyplot as plt
import numpy as np

rf = RandomForestClassifier(n_estimators=100, random_state=42)
rf.fit(X_train, y_train)

probs_uncalibrated = rf.predict_proba(X_test)[:, 1]

calibrated_clf = CalibratedClassifierCV(rf, cv=5, method='sigmoid')
calibrated_clf.fit(X_train, y_train)
probs_calibrated = calibrated_clf.predict_proba(X_test)[:, 1]

frac_pos, mean_pred = calibration_curve(y_test, probs_uncalibrated, n_bins=10)
frac_pos_cal, mean_pred_cal = calibration_curve(y_test, probs_calibrated, n_bins=10)

plt.figure(figsize=(8, 6))
plt.plot(mean_pred, frac_pos, 's-', label='Uncalibrated')
plt.plot(mean_pred_cal, frac_pos_cal, 'o-', label='Calibrated')
plt.plot([0, 1], [0, 1], 'k--', label='Perfect calibration')
plt.xlabel('Mean Predicted Probability')
plt.ylabel('Fraction of Positives')
plt.title('Calibration Curve')
plt.legend()
plt.show()

import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split, cross_val_score, StratifiedKFold, RandomizedSearchCV
from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import StandardScaler, OneHotEncoder
from sklearn.impute import SimpleImputer
from sklearn.ensemble import RandomForestClassifier, GradientBoostingClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import (classification_report, confusion_matrix,
                             roc_auc_score, average_precision_score,
                             precision_recall_curve)
from sklearn.inspection import permutation_importance
from scipy.stats import randint

# 1. Simulate realistic churn data
np.random.seed(42)
n = 5000
df = pd.DataFrame({
    'age': np.random.normal(40, 12, n).clip(18, 80),
    'tenure_months': np.random.exponential(24, n).clip(1, 120),
    'monthly_charges': np.random.normal(65, 20, n).clip(20, 150),
    'total_charges': np.random.normal(4000, 3000, n).clip(0, 20000),
    'num_support_tickets': np.random.poisson(1.5, n),
    'contract_type': np.random.choice(['month-to-month', 'one-year', 'two-year'], n, p=[0.5, 0.3, 0.2]),
    'payment_method': np.random.choice(['credit_card', 'bank_transfer', 'electronic_check', 'mailed_check'], n),
    'internet_service': np.random.choice(['fiber', 'dsl', 'none'], n, p=[0.4, 0.4, 0.2]),
    'has_dependents': np.random.choice([0, 1], n, p=[0.7, 0.3]),
})
for col in ['age', 'monthly_charges', 'num_support_tickets']:
    df.loc[np.random.choice(n, int(n * 0.05), replace=False), col] = np.nan
churn_prob = 0.1 + 0.3 * (df['contract_type'] == 'month-to-month').astype(float)
df['churned'] = (np.random.random(n) < churn_prob).astype(int)

print('Class distribution:', df['churned'].value_counts().to_dict())

# 2. Split
X = df.drop('churned', axis=1)
y = df['churned']
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42, stratify=y)

# 3. Preprocessing
numeric_features = ['age', 'tenure_months', 'monthly_charges', 'total_charges', 'num_support_tickets', 'has_dependents']
categorical_features = ['contract_type', 'payment_method', 'internet_service']

preprocessor = ColumnTransformer([
    ('num', Pipeline([('imputer', SimpleImputer(strategy='median')), ('scaler', StandardScaler())]), numeric_features),
    ('cat', Pipeline([('imputer', SimpleImputer(strategy='most_frequent')), ('onehot', OneHotEncoder(handle_unknown='ignore'))]), categorical_features),
])

# 4. Compare models
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
models = {
    'LogisticRegression': Pipeline([('pre', preprocessor), ('clf', LogisticRegression(max_iter=1000, class_weight='balanced'))]),
    'RandomForest': Pipeline([('pre', preprocessor), ('clf', RandomForestClassifier(200, class_weight='balanced', random_state=42))]),
    'GradientBoosting': Pipeline([('pre', preprocessor), ('clf', GradientBoostingClassifier(200, random_state=42))]),
}

for name, model in models.items():
    f1 = cross_val_score(model, X_train, y_train, cv=cv, scoring='f1')
    auc = cross_val_score(model, X_train, y_train, cv=cv, scoring='roc_auc')
    print(f'{name:20s} F1: {f1.mean():.4f} (+/- {f1.std():.4f})  AUC: {auc.mean():.4f}')

# 5. Tune best model
param_dist = {
    'clf__n_estimators': randint(100, 500),
    'clf__max_depth': [None, 5, 10, 15, 20],
    'clf__max_features': ['sqrt', 'log2'],
    'clf__min_samples_split': randint(2, 20),
}
search = RandomizedSearchCV(models['RandomForest'], param_dist, n_iter=30, cv=cv, scoring='f1', n_jobs=-1, random_state=42)
search.fit(X_train, y_train)
print(f'Best params: {search.best_params_}')

# 6. Evaluate on test set
y_pred = search.predict(X_test)
y_probs = search.predict_proba(X_test)[:, 1]
print(classification_report(y_test, y_pred))
print(f'ROC-AUC: {roc_auc_score(y_test, y_probs):.4f}')
print(f'PR-AUC:  {average_precision_score(y_test, y_probs):.4f}')

# 7. Tune threshold
precisions, recalls, thresholds = precision_recall_curve(y_test, y_probs)
f1_scores = 2 * (precisions * recalls) / (precisions + recalls + 1e-8)
best_threshold = thresholds[np.argmax(f1_scores)]
y_pred_tuned = (y_probs >= best_threshold).astype(int)
print(f'Optimal threshold: {best_threshold:.3f}')
print(classification_report(y_test, y_pred_tuned))

# 8. Feature importance
perm_imp = permutation_importance(search.best_estimator_, X_test, y_test, n_repeats=10, random_state=42, scoring='f1')
feature_names = search.best_estimator_.named_steps['pre'].get_feature_names_out()
imp_df = pd.DataFrame({'feature': feature_names, 'importance': perm_imp.importances_mean}).sort_values('importance', ascending=False)
print(imp_df.head(10).to_string(index=False))

import joblib
import numpy as np

# Save the full pipeline + threshold
model_artifact = {
    'pipeline': search.best_estimator_,
    'threshold': best_threshold,
    'version': '1.0.0',
    'feature_names': list(X.columns),
    'training_date': '2026-03-15',
}
joblib.dump(model_artifact, 'churn_model_v1.joblib', compress=3)

# Load and predict in production
artifact = joblib.load('churn_model_v1.joblib')
pipeline = artifact['pipeline']
threshold = artifact['threshold']

# New data (same schema as training)
new_customer = pd.DataFrame({
    'age': [32],
    'tenure_months': [6],
    'monthly_charges': [85.0],
    'total_charges': [510.0],
    'num_support_tickets': [3],
    'contract_type': ['month-to-month'],
    'payment_method': ['electronic_check'],
    'internet_service': ['fiber'],
    'has_dependents': [0],
})

proba = pipeline.predict_proba(new_customer)[:, 1][0]
prediction = int(proba >= threshold)
print(f'Churn probability: {proba:.3f}')
print(f'Prediction (threshold={threshold:.3f}): {"churn" if prediction else "stay"}')