{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Week 6 Lab: Model Comparisons\n",
    "\n",
    "This lab is about the *tools* for comparing models, in two flavors:\n",
    "\n",
    "- **Part 1 (Regression):** compare four delay-discounting models -- exponential,\n",
    "  hyperbolic, Green & Myerson's hyperboloid, and Rachlin's hyperboloid -- using\n",
    "  R-squared, MAE, RMSE, AIC, BIC, and AICc.\n",
    "- **Part 2 (Classification):** compare logistic regression, a decision tree, a\n",
    "  random forest, and an SVM at predicting challenging behavior, using accuracy,\n",
    "  precision, recall, F1, MCC, and ROC-AUC.\n",
    "\n",
    "Information criteria (AIC/BIC/AICc) matter when models differ in parameter count:\n",
    "they penalize complexity so the comparison is fair."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup\n",
    "\n",
    "Import everything you need: `scipy.optimize` (fitting), `sklearn.metrics`, and the four classifier families from scikit-learn."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import differential_evolution\n",
    "\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.metrics import (accuracy_score, precision_score, recall_score,\n",
    "                             f1_score, matthews_corrcoef, roc_auc_score, roc_curve)\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.tree import DecisionTreeClassifier\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.svm import SVC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 1: Discounting Model Comparison (Regression)\n",
    "\n",
    "## Task 1: Load the discounting data\n",
    "\n",
    "Each row is one participant/commodity with indifference points at 7 delays plus the `Amount` and `Commodity`."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "raw_data = pd.read_csv('participant_discounting_data.csv')\n",
    "raw_data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 2: Define the four discounting models\n",
    "\n",
    "Each returns predicted value as a proportion of the amount (so indifference points are normalized to 0-1).\n",
    "\n",
    "- Exponential: $V = e^{-kD}$\n",
    "- Hyperbolic (Mazur): $V = 1/(1+kD)$\n",
    "- Green-Myerson: $V = 1/(1+kD)^s$\n",
    "- Rachlin: $V = 1/(1+kD^s)$"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "def exponential_model(delay, k):\n",
    "    return np.exp(-k * delay)\n",
    "\n",
    "def hyperbolic_model(delay, k):\n",
    "    return 1 / (1 + k * delay)\n",
    "\n",
    "def green_myerson_hyperboloid(delay, k, s):\n",
    "    return 1 / ((1 + k * delay) ** s)\n",
    "\n",
    "def rachlin_hyperboloid(delay, k, s):\n",
    "    return 1 / (1 + k * (delay ** s))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 3: Helpers to fit a model and compute fit metrics\n",
    "\n",
    "`fit_model_robust` uses global optimization (`differential_evolution`) to minimize squared error. `calculate_fit_metrics` returns R-squared, MAE, RMSE plus the information criteria AIC, BIC, and AICc (small-sample corrected)."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "def prepare_data_for_fitting(raw_data, participant_id, commodity):\n",
    "    row = raw_data[(raw_data['participant_id'] == participant_id) &\n",
    "                   (raw_data['Commodity'] == commodity)].iloc[0]\n",
    "    delay_cols = [c for c in raw_data.columns if c.startswith('delay_')]\n",
    "    delays = np.array([float(c.split('_')[1]) for c in delay_cols])\n",
    "    indiff = row[delay_cols].values.astype(float) / row['Amount']   # normalize to 0-1\n",
    "    return delays, indiff\n",
    "\n",
    "def fit_model_robust(model_func, delays, indiff, bounds):\n",
    "    def objective(params):\n",
    "        try:\n",
    "            return np.sum((indiff - model_func(delays, *params)) ** 2)\n",
    "        except Exception:\n",
    "            return 1e10\n",
    "    result = differential_evolution(objective, bounds, seed=0, maxiter=200, tol=1e-6)\n",
    "    params = result.x\n",
    "    predictions = model_func(delays, *params)\n",
    "    return params, predictions\n",
    "\n",
    "def calculate_fit_metrics(observed, predicted, n_params, n_obs):\n",
    "    residuals = observed - predicted\n",
    "    ss_res = np.sum(residuals ** 2)\n",
    "    ss_tot = np.sum((observed - np.mean(observed)) ** 2)\n",
    "    r_squared = 1 - (ss_res / ss_tot) if ss_tot > 0 else 0\n",
    "    mae = np.mean(np.abs(residuals))\n",
    "    rmse = np.sqrt(np.mean(residuals ** 2))\n",
    "    mse = ss_res / n_obs\n",
    "    log_lik = -0.5 * n_obs * (np.log(2 * np.pi * mse) + 1) if mse > 0 else 0\n",
    "    aic = 2 * n_params - 2 * log_lik\n",
    "    bic = np.log(n_obs) * n_params - 2 * log_lik\n",
    "    aicc = aic + (2 * n_params * (n_params + 1)) / (n_obs - n_params - 1) \\\n",
    "        if n_obs - n_params - 1 > 0 else np.inf\n",
    "    return {'R2': r_squared, 'MAE': mae, 'RMSE': rmse, 'AIC': aic, 'BIC': bic, 'AICc': aicc}\n",
    "\n",
    "MODELS = {\n",
    "    'Exponential':   {'func': exponential_model,        'bounds': [(1e-6, 10)],            'n_params': 1},\n",
    "    'Hyperbolic':    {'func': hyperbolic_model,         'bounds': [(1e-6, 10)],            'n_params': 1},\n",
    "    'Green-Myerson': {'func': green_myerson_hyperboloid,'bounds': [(1e-6, 10), (0.1, 5)],  'n_params': 2},\n",
    "    'Rachlin':       {'func': rachlin_hyperboloid,      'bounds': [(1e-6, 10), (0.1, 5)],  'n_params': 2},\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 4: Fit every model to every participant\n",
    "\n",
    "Loop over participants, fit all four models, and collect the metrics into a tidy summary DataFrame. (Fitting all participants on 7 points each takes a moment.)"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "summary_rows = []\n",
    "for pid in raw_data['participant_id'].unique():\n",
    "    for comm in raw_data[raw_data['participant_id'] == pid]['Commodity'].unique():\n",
    "        delays, indiff = prepare_data_for_fitting(raw_data, pid, comm)\n",
    "        for name, m in MODELS.items():\n",
    "            params, preds = fit_model_robust(m['func'], delays, indiff, m['bounds'])\n",
    "            metrics = calculate_fit_metrics(indiff, preds, m['n_params'], len(indiff))\n",
    "            summary_rows.append({'Participant': pid, 'Commodity': comm, 'Model': name, **metrics})\n",
    "\n",
    "summary_df = pd.DataFrame(summary_rows)\n",
    "summary_df.head(8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 5: Compare the models across metrics\n",
    "\n",
    "For each fit metric, show its distribution by model so you can see which model wins on which metric. Note how the information criteria treat the 1- vs 2-parameter models differently."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "mean_by_model = summary_df.groupby('Model')[['R2', 'MAE', 'RMSE', 'AIC', 'BIC', 'AICc']].mean()\n",
    "print(mean_by_model.round(4))\n",
    "\n",
    "metrics_to_plot = ['R2', 'MAE', 'RMSE', 'AIC', 'BIC', 'AICc']\n",
    "fig, axes = plt.subplots(2, 3, figsize=(16, 9))\n",
    "for ax, metric in zip(axes.ravel(), metrics_to_plot):\n",
    "    sns.boxplot(data=summary_df, x='Model', y=metric, ax=ax)\n",
    "    ax.set_title(metric)\n",
    "    ax.tick_params(axis='x', rotation=30)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 2: Challenging Behavior Prediction (Classification)\n",
    "\n",
    "## Task 6: Load and prepare the data\n",
    "\n",
    "Load `challenging_behavior_data.csv`, split features/label, make a train/test split, and standardize features (needed by logistic regression and the SVM). The SVM is expensive on large data, so we sample a subset for training speed."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "clf_data = pd.read_csv('challenging_behavior_data.csv')\n",
    "\n",
    "# RBF-SVM training cost grows steeply with n; sample for a tractable demo.\n",
    "clf_data = clf_data.sample(n=5000, random_state=42).reset_index(drop=True)\n",
    "\n",
    "X = clf_data.drop('ChallengingBehavior', axis=1)\n",
    "y = clf_data['ChallengingBehavior']\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)\n",
    "\n",
    "scaler = StandardScaler()\n",
    "X_train_scaled = scaler.fit_transform(X_train)\n",
    "X_test_scaled = scaler.transform(X_test)\n",
    "print('train:', X_train.shape, ' test:', X_test.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 7: Fit four classifiers and score them\n",
    "\n",
    "Fit logistic regression, a decision tree, a random forest, and an RBF SVM. For each, compute accuracy, precision, recall, F1, MCC, and ROC-AUC, and add its ROC curve to a shared plot. Scale-sensitive models (logistic, SVM) use the scaled features."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "classifiers = {\n",
    "    'Logistic Regression': LogisticRegression(max_iter=1000),\n",
    "    'Decision Tree': DecisionTreeClassifier(random_state=42),\n",
    "    'Random Forest': RandomForestClassifier(random_state=42),\n",
    "    'SVM (RBF)': SVC(probability=True, random_state=42),\n",
    "}\n",
    "\n",
    "results = []\n",
    "plt.figure(figsize=(9, 7))\n",
    "for name, model in classifiers.items():\n",
    "    if 'SVM' in name or 'Logistic' in name:\n",
    "        model.fit(X_train_scaled, y_train)\n",
    "        y_pred = model.predict(X_test_scaled)\n",
    "        y_proba = model.predict_proba(X_test_scaled)[:, 1]\n",
    "    else:\n",
    "        model.fit(X_train, y_train)\n",
    "        y_pred = model.predict(X_test)\n",
    "        y_proba = model.predict_proba(X_test)[:, 1]\n",
    "\n",
    "    results.append({\n",
    "        'Model': name,\n",
    "        'Accuracy': accuracy_score(y_test, y_pred),\n",
    "        'Precision': precision_score(y_test, y_pred),\n",
    "        'Recall': recall_score(y_test, y_pred),\n",
    "        'F1 Score': f1_score(y_test, y_pred),\n",
    "        'MCC': matthews_corrcoef(y_test, y_pred),\n",
    "        'ROC AUC': roc_auc_score(y_test, y_proba),\n",
    "    })\n",
    "    fpr, tpr, _ = roc_curve(y_test, y_proba)\n",
    "    plt.plot(fpr, tpr, label=f\"{name} (AUC = {results[-1]['ROC AUC']:.2f})\")\n",
    "\n",
    "plt.plot([0, 1], [0, 1], 'k--', label='Random Guess')\n",
    "plt.xlabel('False Positive Rate'); plt.ylabel('True Positive Rate')\n",
    "plt.title('ROC Curves'); plt.legend(); plt.grid(True); plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 8: Results table and interpretation\n",
    "\n",
    "Put the classification metrics in one table. Which model would you choose, and on which metric -- and why might accuracy alone be misleading if challenging behavior is rare?"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "results_df = pd.DataFrame(results).set_index('Model').round(3)\n",
    "print(results_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Wrap-up\n",
    "\n",
    "Across both parts: when does a more complex model (extra parameter, ensemble) actually earn its complexity? Tie your answer to what the information criteria showed in Part 1 and to the precision/recall trade-off in Part 2."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}