{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.11","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":102335,"databundleVersionId":12518947,"sourceType":"competition"}],"dockerImageVersionId":31040,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# **CMI BFRB Behavior Detection - Comprehensive EDA**\n\n## Overview\nThis is an in-depth exploratory data analysis (EDA) on the CMI BFRB (Body-Focused Repetitive Behavior) detection dataset, focusing on sensor-based gesture recognition and behavior classification.\n\n## Dataset Description\nThe dataset contains multimodal sensor data from wearable devices designed to detect repetitive behaviors like hair pulling, skin picking, and nail biting. Key components include:\n\n- **Sensor Data**: IMU (accelerometer/gyroscope), thermopile, and time-of-flight sensors\n- **Behavioral Phases**: Each sequence contains Transition → Pause → Gesture phases\n- **Demographics**: Subject characteristics including age, height, handedness, and body measurements\n- **Gesture Types**: Target (BFRB behaviors) and Non-target sequences\n\n## Analysis Components\n\n### 1. **Data Structure Analysis**\n- Missing data patterns across sensor modalities\n- Sequence length distributions and characteristics\n- Subject-level statistics and demographics correlation\n\n### 2. **Sensor-Specific Analysis**\n- **IMU Data**: Acceleration and rotation patterns with behavior phase markers\n- **Thermopile**: Temperature sensor availability and correlation patterns\n- **Time-of-Flight**: Spatial distance measurements and grid analysis\n\n### 3. **Behavioral Pattern Analysis**\n- Phase duration statistics and distributions\n- Gesture-specific sensor signatures\n- Subject variability and demographic influences\n\n### 4. **Advanced Feature Engineering**\n- **Phase-specific features**: Separate statistics for each behavioral phase\n- **Motion dynamics**: Velocity, jerk, and energy calculations\n- **Frequency domain**: FFT-based spectral analysis\n- **Spatial features**: Center of mass and variance for ToF sensor grids\n- **Statistical discriminators**: ANOVA-based feature ranking\n\n### 5. **Signal Processing**\n- Time and frequency domain analysis\n- Moving statistics and behavior phase transitions\n- Correlation analysis across sensor modalities\n\n## Key Insights\n- Gesture-phase features show strongest discriminative power\n- Motion magnitude features outperform individual axis measurements\n- Phase-specific analysis reveals distinct behavioral signatures\n- Multi-sensor fusion provides complementary information for classification\n\n## Applications\nThis analysis provides the foundation for building machine learning models to automatically detect BFRB behaviors in real-time, enabling therapeutic interventions and behavioral monitoring.","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom scipy import stats\nfrom scipy.fft import fft, fftfreq\nimport plotly.express as px\nimport plotly.graph_objects as go\nfrom plotly.subplots import make_subplots\nimport warnings\nwarnings.filterwarnings('ignore')\n\n# Set up plotting style\nplt.style.use('default')\nsns.set_palette(\"husl\")\nplt.rcParams['figure.figsize'] = (12, 8)\n\n# Load data\nprint(\"Loading data...\")\ntrain_df = pd.read_csv(\"/kaggle/input/cmi-detect-behavior-with-sensor-data/train.csv\")\ntest_df = pd.read_csv(\"/kaggle/input/cmi-detect-behavior-with-sensor-data/test.csv\")\ntrain_demographics_df = pd.read_csv(\"/kaggle/input/cmi-detect-behavior-with-sensor-data/train_demographics.csv\")\ntest_demographics_df = pd.read_csv(\"/kaggle/input/cmi-detect-behavior-with-sensor-data/test_demographics.csv\")\n\nprint(f\"Train shape: {train_df.shape}\")\nprint(f\"Test shape: {test_df.shape}\")\nprint(f\"Train demographics shape: {train_demographics_df.shape}\")\nprint(f\"Test demographics shape: {test_demographics_df.shape}\")\n\n# ================================\n# 1. BASIC DATA OVERVIEW\n# ================================\n\nprint(\"\\n\" + \"=\"*50)\nprint(\"1. BASIC DATA OVERVIEW\")\nprint(\"=\"*50)\n\n# Basic info\nprint(\"\\nTrain DataFrame Info:\")\nprint(train_df.info())\n\nprint(\"\\nUnique sequence types:\")\nprint(train_df['sequence_type'].value_counts())\n\nprint(\"\\nUnique behaviors:\")\nprint(train_df['behavior'].value_counts())\n\nprint(\"\\nUnique gestures:\")\nprint(train_df['gesture'].value_counts())\n\n# Missing data analysis\nprint(\"\\nMissing data analysis:\")\nmissing_data = train_df.isnull().sum().sort_values(ascending=False)\nmissing_percentage = (missing_data / len(train_df)) * 100\nmissing_df = pd.DataFrame({\n    'Missing_Count': missing_data,\n    'Percentage': missing_percentage\n})\nprint(missing_df[missing_df['Missing_Count'] > 0].head(20))\n\n# ================================\n# 2. SEQUENCE-LEVEL ANALYSIS\n# ================================\n\nprint(\"\\n\" + \"=\"*50)\nprint(\"2. SEQUENCE-LEVEL ANALYSIS\")\nprint(\"=\"*50)\n\n# Sequence statistics\nsequence_stats = train_df.groupby('sequence_id').agg({\n    'sequence_counter': 'count',\n    'subject': 'first',\n    'gesture': 'first',\n    'sequence_type': 'first',\n    'behavior': lambda x: x.value_counts().to_dict()\n}).reset_index()\n\nsequence_stats.columns = ['sequence_id', 'sequence_length', 'subject', 'gesture', 'sequence_type', 'behavior_counts']\n\nprint(f\"Total sequences: {len(sequence_stats)}\")\nprint(f\"Sequence length statistics:\")\nprint(sequence_stats['sequence_length'].describe())\n\n# Plot sequence length distribution\nfig, axes = plt.subplots(2, 2, figsize=(15, 12))\n\n# Sequence length distribution\naxes[0,0].hist(sequence_stats['sequence_length'], bins=50, alpha=0.7, edgecolor='black')\naxes[0,0].set_title('Distribution of Sequence Lengths')\naxes[0,0].set_xlabel('Sequence Length')\naxes[0,0].set_ylabel('Frequency')\n\n# Sequence length by gesture type\ntarget_sequences = sequence_stats[sequence_stats['sequence_type'] == 'Target']\nnon_target_sequences = sequence_stats[sequence_stats['sequence_type'] == 'Non-target']\n\naxes[0,1].hist([target_sequences['sequence_length'], non_target_sequences['sequence_length']], \n               bins=30, alpha=0.7, label=['Target (BFRB)', 'Non-target'], edgecolor='black')\naxes[0,1].set_title('Sequence Length by Type')\naxes[0,1].set_xlabel('Sequence Length')\naxes[0,1].set_ylabel('Frequency')\naxes[0,1].legend()\n\n# Subjects distribution\naxes[1,0].hist(sequence_stats.groupby('subject').size(), bins=20, alpha=0.7, edgecolor='black')\naxes[1,0].set_title('Sequences per Subject')\naxes[1,0].set_xlabel('Number of Sequences')\naxes[1,0].set_ylabel('Number of Subjects')\n\n# Gesture distribution\ngesture_counts = sequence_stats['gesture'].value_counts()\ntop_gestures = gesture_counts.head(10)\naxes[1,1].barh(range(len(top_gestures)), top_gestures.values)\naxes[1,1].set_yticks(range(len(top_gestures)))\naxes[1,1].set_yticklabels(top_gestures.index, fontsize=8)\naxes[1,1].set_title('Top 10 Gestures Distribution')\naxes[1,1].set_xlabel('Count')\n\nplt.tight_layout()\nplt.show()\n\n# ================================\n# 3. DEMOGRAPHICS ANALYSIS\n# ================================\n\nprint(\"\\n\" + \"=\"*50)\nprint(\"3. DEMOGRAPHICS ANALYSIS\")\nprint(\"=\"*50)\n\nprint(\"Demographics summary:\")\nprint(train_demographics_df.describe())\n\n# Plot demographics\nfig, axes = plt.subplots(2, 3, figsize=(18, 12))\n\n# Age distribution\naxes[0,0].hist(train_demographics_df['age'], bins=20, alpha=0.7, edgecolor='black')\naxes[0,0].set_title('Age Distribution')\naxes[0,0].set_xlabel('Age')\naxes[0,0].set_ylabel('Frequency')\n\n# Adult vs Child\nadult_child_counts = train_demographics_df['adult_child'].value_counts()\naxes[0,1].bar(['Child', 'Adult'], adult_child_counts.values, alpha=0.7, edgecolor='black')\naxes[0,1].set_title('Adult vs Child Distribution')\naxes[0,1].set_ylabel('Count')\n\n# Sex distribution\nsex_counts = train_demographics_df['sex'].value_counts()\naxes[0,2].bar(['Female', 'Male'], sex_counts.values, alpha=0.7, edgecolor='black')\naxes[0,2].set_title('Sex Distribution')\naxes[0,2].set_ylabel('Count')\n\n# Height distribution\naxes[1,0].hist(train_demographics_df['height_cm'], bins=20, alpha=0.7, edgecolor='black')\naxes[1,0].set_title('Height Distribution')\naxes[1,0].set_xlabel('Height (cm)')\naxes[1,0].set_ylabel('Frequency')\n\n# Arm measurements\naxes[1,1].scatter(train_demographics_df['shoulder_to_wrist_cm'], \n                  train_demographics_df['elbow_to_wrist_cm'], alpha=0.6)\naxes[1,1].set_title('Arm Measurements Correlation')\naxes[1,1].set_xlabel('Shoulder to Wrist (cm)')\naxes[1,1].set_ylabel('Elbow to Wrist (cm)')\n\n# Handedness\nhandedness_counts = train_demographics_df['handedness'].value_counts()\naxes[1,2].bar(['Left', 'Right'], handedness_counts.values, alpha=0.7, edgecolor='black')\naxes[1,2].set_title('Handedness Distribution')\naxes[1,2].set_ylabel('Count')\n\nplt.tight_layout()\nplt.show()\n\n# ================================\n# 4. SENSOR DATA ANALYSIS\n# ================================\n\nprint(\"\\n\" + \"=\"*50)\nprint(\"4. SENSOR DATA ANALYSIS\")\nprint(\"=\"*50)\n\n# Define sensor columns\nimu_cols = ['acc_x', 'acc_y', 'acc_z', 'rot_w', 'rot_x', 'rot_y', 'rot_z']\nthm_cols = [col for col in train_df.columns if 'thm_' in col]\ntof_cols = [col for col in train_df.columns if 'tof_' in col]\n\nprint(f\"IMU columns: {len(imu_cols)}\")\nprint(f\"Thermopile columns: {len(thm_cols)}\")\nprint(f\"Time-of-flight columns: {len(tof_cols)}\")\n\n# Sensor availability analysis\nprint(\"\\nSensor availability analysis:\")\nfor sensor_group, cols in [('IMU', imu_cols), ('Thermopile', thm_cols), ('ToF', tof_cols)]:\n    if cols:\n        available_data = train_df[cols].notna().all(axis=1).sum()\n        total_rows = len(train_df)\n        print(f\"{sensor_group}: {available_data}/{total_rows} ({available_data/total_rows*100:.1f}%) complete rows\")\n\n# Sample a few sequences for detailed analysis\nsample_sequences = train_df['sequence_id'].unique()[:5]\nprint(f\"\\nAnalyzing sample sequences: {sample_sequences}\")\n\n# Plot IMU data for sample sequences\nfig, axes = plt.subplots(len(sample_sequences), 3, figsize=(20, 4*len(sample_sequences)))\nif len(sample_sequences) == 1:\n    axes = axes.reshape(1, -1)\n\nfor i, seq_id in enumerate(sample_sequences):\n    seq_data = train_df[train_df['sequence_id'] == seq_id].reset_index(drop=True)\n    gesture = seq_data['gesture'].iloc[0]\n    \n    # Acceleration\n    axes[i,0].plot(seq_data['acc_x'], label='X', alpha=0.7)\n    axes[i,0].plot(seq_data['acc_y'], label='Y', alpha=0.7)\n    axes[i,0].plot(seq_data['acc_z'], label='Z', alpha=0.7)\n    axes[i,0].set_title(f'Acceleration - {gesture}')\n    axes[i,0].set_ylabel('Acceleration (m/s²)')\n    axes[i,0].legend()\n    axes[i,0].grid(True, alpha=0.3)\n    \n    # Add behavior phase markers\n    behavior_changes = seq_data[seq_data['behavior'] != seq_data['behavior'].shift()].index\n    for change_idx in behavior_changes:\n        if change_idx > 0:\n            axes[i,0].axvline(x=change_idx, color='red', linestyle='--', alpha=0.5)\n    \n    # Rotation\n    axes[i,1].plot(seq_data['rot_w'], label='W', alpha=0.7)\n    axes[i,1].plot(seq_data['rot_x'], label='X', alpha=0.7)\n    axes[i,1].plot(seq_data['rot_y'], label='Y', alpha=0.7)\n    axes[i,1].plot(seq_data['rot_z'], label='Z', alpha=0.7)\n    axes[i,1].set_title(f'Rotation - {gesture}')\n    axes[i,1].set_ylabel('Rotation (quaternion)')\n    axes[i,1].legend()\n    axes[i,1].grid(True, alpha=0.3)\n    \n    # Add behavior phase markers\n    for change_idx in behavior_changes:\n        if change_idx > 0:\n            axes[i,1].axvline(x=change_idx, color='red', linestyle='--', alpha=0.5)\n    \n    # Acceleration magnitude\n    acc_magnitude = np.sqrt(seq_data['acc_x']**2 + seq_data['acc_y']**2 + seq_data['acc_z']**2)\n    axes[i,2].plot(acc_magnitude, color='purple', alpha=0.7)\n    axes[i,2].set_title(f'Acceleration Magnitude - {gesture}')\n    axes[i,2].set_ylabel('Magnitude (m/s²)')\n    axes[i,2].grid(True, alpha=0.3)\n    \n    # Add behavior phase markers\n    for change_idx in behavior_changes:\n        if change_idx > 0:\n            axes[i,2].axvline(x=change_idx, color='red', linestyle='--', alpha=0.5)\n\nplt.tight_layout()\nplt.show()\n\n# ================================\n# 5. BEHAVIOR PHASE ANALYSIS\n# ================================\n\nprint(\"\\n\" + \"=\"*50)\nprint(\"5. BEHAVIOR PHASE ANALYSIS\")\nprint(\"=\"*50)\n\n# Analyze behavior phases\nphase_stats = train_df.groupby(['sequence_id', 'behavior']).size().reset_index(name='count')\nphase_summary = phase_stats.groupby('behavior')['count'].agg(['mean', 'std', 'min', 'max']).round(2)\nprint(\"Behavior phase duration statistics:\")\nprint(phase_summary)\n\n# Plot behavior phase durations\nplt.figure(figsize=(15, 5))\n\nplt.subplot(1, 3, 1)\nphase_counts = train_df['behavior'].value_counts()\nplt.bar(phase_counts.index, phase_counts.values, alpha=0.7, edgecolor='black')\nplt.title('Behavior Phase Distribution')\nplt.ylabel('Count')\nplt.xticks(rotation=45)\n\nplt.subplot(1, 3, 2)\nsns.boxplot(data=phase_stats, x='behavior', y='count')\nplt.title('Behavior Phase Duration Distribution')\nplt.xticks(rotation=45)\n\nplt.subplot(1, 3, 3)\nfor behavior in phase_stats['behavior'].unique():\n    behavior_data = phase_stats[phase_stats['behavior'] == behavior]['count']\n    plt.hist(behavior_data, alpha=0.6, label=behavior, bins=20)\nplt.title('Behavior Phase Duration Histograms')\nplt.xlabel('Duration (data points)')\nplt.ylabel('Frequency')\nplt.legend()\n\nplt.tight_layout()\nplt.show()\n\n# ================================\n# 6. GESTURE-SPECIFIC ANALYSIS\n# ================================\n\nprint(\"\\n\" + \"=\"*50)\nprint(\"6. GESTURE-SPECIFIC ANALYSIS\")\nprint(\"=\"*50)\n\n# Focus on target sequences (BFRB gestures)\ntarget_data = train_df[train_df['sequence_type'] == 'Target'].copy()\nprint(f\"Target sequences: {len(target_data['sequence_id'].unique())}\")\n\n# Gesture statistics\ngesture_stats = target_data.groupby('gesture').agg({\n    'sequence_id': 'nunique',\n    'subject': 'nunique'\n}).reset_index()\ngesture_stats.columns = ['gesture', 'num_sequences', 'num_subjects']\n\nprint(\"\\nGesture statistics:\")\nprint(gesture_stats.sort_values('num_sequences', ascending=False))\n\n# Plot gesture distribution\nplt.figure(figsize=(15, 8))\ngesture_counts = target_data['gesture'].value_counts()\nplt.barh(range(len(gesture_counts)), gesture_counts.values)\nplt.yticks(range(len(gesture_counts)), gesture_counts.index)\nplt.title('BFRB Gesture Distribution')\nplt.xlabel('Number of Data Points')\nplt.tight_layout()\nplt.show()\n\n# ================================\n# 7. THERMOPILE ANALYSIS\n# ================================\n\nprint(\"\\n\" + \"=\"*50)\nprint(\"7. THERMOPILE ANALYSIS\")\nprint(\"=\"*50)\n\nif thm_cols:\n    # Thermopile availability\n    thm_available = target_data[thm_cols].notna().any(axis=1)\n    print(f\"Sequences with thermopile data: {thm_available.sum()}/{len(target_data)} ({thm_available.sum()/len(target_data)*100:.1f}%)\")\n    \n    # Thermopile statistics\n    thm_data = target_data[thm_cols].dropna()\n    if not thm_data.empty:\n        print(\"\\nThermopile statistics:\")\n        print(thm_data.describe())\n        \n        # Plot thermopile distributions\n        fig, axes = plt.subplots(2, 3, figsize=(18, 12))\n        axes = axes.flatten()\n        \n        for i, col in enumerate(thm_cols):\n            if i < len(axes):\n                thm_data[col].hist(bins=50, alpha=0.7, ax=axes[i])\n                axes[i].set_title(f'{col} Distribution')\n                axes[i].set_xlabel('Temperature (°C)')\n                axes[i].set_ylabel('Frequency')\n        \n        # Remove empty subplots\n        for i in range(len(thm_cols), len(axes)):\n            fig.delaxes(axes[i])\n        \n        plt.tight_layout()\n        plt.show()\n        \n        # Thermopile correlation matrix\n        plt.figure(figsize=(8, 6))\n        thm_corr = thm_data.corr()\n        sns.heatmap(thm_corr, annot=True, cmap='coolwarm', center=0)\n        plt.title('Thermopile Sensors Correlation Matrix')\n        plt.tight_layout()\n        plt.show()\n\n# ================================\n# 8. TIME-OF-FLIGHT ANALYSIS\n# ================================\n\nprint(\"\\n\" + \"=\"*50)\nprint(\"8. TIME-OF-FLIGHT ANALYSIS\")\nprint(\"=\"*50)\n\nif tof_cols:\n    # ToF availability (excluding -1 values which indicate no response)\n    tof_data = target_data[tof_cols].replace(-1, np.nan)\n    tof_available = tof_data.notna().any(axis=1)\n    print(f\"Sequences with ToF data: {tof_available.sum()}/{len(target_data)} ({tof_available.sum()/len(target_data)*100:.1f}%)\")\n    \n    # ToF statistics (excluding -1 values)\n    tof_valid = tof_data.dropna()\n    if not tof_valid.empty:\n        print(\"\\nToF statistics (valid readings only):\")\n        print(tof_valid.describe())\n        \n        # Analyze ToF sensor patterns\n        # Group ToF columns by sensor\n        tof_sensors = {}\n        for col in tof_cols:\n            sensor_num = col.split('_')[1]\n            if sensor_num not in tof_sensors:\n                tof_sensors[sensor_num] = []\n            tof_sensors[sensor_num].append(col)\n        \n        print(f\"\\nToF sensors found: {list(tof_sensors.keys())}\")\n        \n        # Plot ToF sensor availability\n        plt.figure(figsize=(12, 6))\n        sensor_availability = []\n        sensor_labels = []\n        \n        for sensor_num, sensor_cols in tof_sensors.items():\n            sensor_data = target_data[sensor_cols].replace(-1, np.nan)\n            availability = sensor_data.notna().any(axis=1).sum() / len(target_data) * 100\n            sensor_availability.append(availability)\n            sensor_labels.append(f'ToF Sensor {sensor_num}')\n        \n        plt.bar(sensor_labels, sensor_availability)\n        plt.title('ToF Sensor Data Availability')\n        plt.ylabel('Availability (%)')\n        plt.xticks(rotation=45)\n        plt.tight_layout()\n        plt.show()\n\n# ================================\n# 9. SUBJECT-SPECIFIC PATTERNS\n# ================================\n\nprint(\"\\n\" + \"=\"*50)\nprint(\"9. SUBJECT-SPECIFIC PATTERNS\")\nprint(\"=\"*50)\n\n# Subject statistics\nsubject_stats = target_data.groupby('subject').agg({\n    'sequence_id': 'nunique',\n    'gesture': lambda x: x.nunique()\n}).reset_index()\nsubject_stats.columns = ['subject', 'num_sequences', 'num_unique_gestures']\n\nprint(\"Subject statistics:\")\nprint(subject_stats.describe())\n\n# Merge with demographics\nsubject_demo = subject_stats.merge(train_demographics_df, on='subject', how='left')\n\n# Plot subject patterns\nfig, axes = plt.subplots(2, 2, figsize=(15, 12))\n\n# Sequences per subject\naxes[0,0].hist(subject_stats['num_sequences'], bins=20, alpha=0.7, edgecolor='black')\naxes[0,0].set_title('Sequences per Subject')\naxes[0,0].set_xlabel('Number of Sequences')\naxes[0,0].set_ylabel('Number of Subjects')\n\n# Unique gestures per subject\naxes[0,1].hist(subject_stats['num_unique_gestures'], bins=10, alpha=0.7, edgecolor='black')\naxes[0,1].set_title('Unique Gestures per Subject')\naxes[0,1].set_xlabel('Number of Unique Gestures')\naxes[0,1].set_ylabel('Number of Subjects')\n\n# Age vs sequences\nif 'age' in subject_demo.columns:\n    axes[1,0].scatter(subject_demo['age'], subject_demo['num_sequences'], alpha=0.6)\n    axes[1,0].set_title('Age vs Number of Sequences')\n    axes[1,0].set_xlabel('Age')\n    axes[1,0].set_ylabel('Number of Sequences')\n\n# Height vs sequences\nif 'height_cm' in subject_demo.columns:\n    axes[1,1].scatter(subject_demo['height_cm'], subject_demo['num_sequences'], alpha=0.6)\n    axes[1,1].set_title('Height vs Number of Sequences')\n    axes[1,1].set_xlabel('Height (cm)')\n    axes[1,1].set_ylabel('Number of Sequences')\n\nplt.tight_layout()\nplt.show()\n\n# ================================\n# 10. SIGNAL PROCESSING ANALYSIS\n# ================================\n\nprint(\"\\n\" + \"=\"*50)\nprint(\"10. SIGNAL PROCESSING ANALYSIS\")\nprint(\"=\"*50)\n\n# Select a representative sequence for signal analysis\nsample_seq_id = target_data['sequence_id'].value_counts().index[0]  # Most common sequence\nsample_data = target_data[target_data['sequence_id'] == sample_seq_id].reset_index(drop=True)\n\nprint(f\"Analyzing sequence {sample_seq_id}: {sample_data['gesture'].iloc[0]}\")\n\n# Calculate acceleration magnitude\nacc_magnitude = np.sqrt(sample_data['acc_x']**2 + sample_data['acc_y']**2 + sample_data['acc_z']**2)\n\n# FFT analysis\nif len(acc_magnitude) > 1:\n    fft_values = np.abs(fft(acc_magnitude))\n    freqs = fftfreq(len(acc_magnitude), d=1.0)  # Assuming 1 Hz sampling for demo\n    \n    plt.figure(figsize=(15, 10))\n    \n    # Time domain\n    plt.subplot(2, 3, 1)\n    plt.plot(acc_magnitude)\n    plt.title('Acceleration Magnitude (Time Domain)')\n    plt.xlabel('Sample')\n    plt.ylabel('Magnitude')\n    \n    # Frequency domain\n    plt.subplot(2, 3, 2)\n    plt.plot(freqs[:len(freqs)//2], fft_values[:len(fft_values)//2])\n    plt.title('Acceleration Magnitude (Frequency Domain)')\n    plt.xlabel('Frequency')\n    plt.ylabel('Magnitude')\n    \n    # Individual acceleration components\n    plt.subplot(2, 3, 3)\n    plt.plot(sample_data['acc_x'], label='X', alpha=0.7)\n    plt.plot(sample_data['acc_y'], label='Y', alpha=0.7)\n    plt.plot(sample_data['acc_z'], label='Z', alpha=0.7)\n    plt.title('Acceleration Components')\n    plt.legend()\n    \n    # Rotation components\n    plt.subplot(2, 3, 4)\n    plt.plot(sample_data['rot_w'], label='W', alpha=0.7)\n    plt.plot(sample_data['rot_x'], label='X', alpha=0.7)\n    plt.plot(sample_data['rot_y'], label='Y', alpha=0.7)\n    plt.plot(sample_data['rot_z'], label='Z', alpha=0.7)\n    plt.title('Rotation Components')\n    plt.legend()\n    \n    # Moving statistics\n    window_size = min(10, len(acc_magnitude)//4)\n    if window_size > 1:\n        moving_mean = pd.Series(acc_magnitude).rolling(window=window_size).mean()\n        moving_std = pd.Series(acc_magnitude).rolling(window=window_size).std()\n        \n        plt.subplot(2, 3, 5)\n        plt.plot(moving_mean, label='Moving Mean', alpha=0.8)\n        plt.plot(moving_std, label='Moving Std', alpha=0.8)\n        plt.title(f'Moving Statistics (window={window_size})')\n        plt.legend()\n    \n    # Behavior phases\n    plt.subplot(2, 3, 6)\n    behavior_numeric = pd.Categorical(sample_data['behavior']).codes\n    plt.plot(behavior_numeric, marker='o')\n    plt.title('Behavior Phases')\n    plt.ylabel('Phase Code')\n    unique_behaviors = sample_data['behavior'].unique()\n    plt.yticks(range(len(unique_behaviors)), unique_behaviors)\n    \n    plt.tight_layout()\n    plt.show()\n\n# ================================\n# 11. CORRELATION ANALYSIS\n# ================================\n\nprint(\"\\n\" + \"=\"*50)\nprint(\"11. CORRELATION ANALYSIS\")\nprint(\"=\"*50)\n\n# Create correlation matrix for IMU data\nimu_data = target_data[imu_cols].corr()\n\nplt.figure(figsize=(10, 8))\nsns.heatmap(imu_data, annot=True, cmap='coolwarm', center=0, fmt='.2f')\nplt.title('IMU Sensors Correlation Matrix')\nplt.tight_layout()\nplt.show()\n\n# ================================\n# 12. DATA QUALITY ASSESSMENT\n# ================================\n\nprint(\"\\n\" + \"=\"*50)\nprint(\"12. DATA QUALITY ASSESSMENT\")\nprint(\"=\"*50)\n\n# Check for anomalies and data quality issues\nprint(\"Data quality assessment:\")\n\n# Check for constant values\nconstant_cols = []\nfor col in imu_cols:\n    if target_data[col].nunique() == 1:\n        constant_cols.append(col)\n\nif constant_cols:\n    print(f\"Columns with constant values: {constant_cols}\")\nelse:\n    print(\"No columns with constant values found\")\n\n# Check for extreme values\nprint(\"\\nExtreme values detection:\")\nfor col in imu_cols:\n    Q1 = target_data[col].quantile(0.25)\n    Q3 = target_data[col].quantile(0.75)\n    IQR = Q3 - Q1\n    lower_bound = Q1 - 1.5 * IQR\n    upper_bound = Q3 + 1.5 * IQR\n    \n    outliers = target_data[(target_data[col] < lower_bound) | (target_data[col] > upper_bound)]\n    if len(outliers) > 0:\n        print(f\"{col}: {len(outliers)} outliers ({len(outliers)/len(target_data)*100:.2f}%)\")\n\n# ================================\n# 13. FEATURE ENGINEERING ANALYSIS\n# ================================\n\nprint(\"\\n\" + \"=\"*50)\nprint(\"13. FEATURE ENGINEERING ANALYSIS\")\nprint(\"=\"*50)\n\ndef extract_sequence_features(df, sequence_id):\n    \"\"\"Extract comprehensive features from a single sequence\"\"\"\n    seq_data = df[df['sequence_id'] == sequence_id].reset_index(drop=True)\n    features = {}\n    \n    # Basic sequence info\n    features['sequence_length'] = len(seq_data)\n    features['gesture'] = seq_data['gesture'].iloc[0]\n    \n    # Phase-specific features\n    for phase in ['Transition', 'Pause', 'Gesture']:\n        phase_data = seq_data[seq_data['behavior'] == phase]\n        if not phase_data.empty:\n            features[f'{phase.lower()}_length'] = len(phase_data)\n            features[f'{phase.lower()}_ratio'] = len(phase_data) / len(seq_data)\n            \n            # IMU features for each phase\n            if all(col in phase_data.columns for col in imu_cols):\n                # Acceleration magnitude\n                acc_mag = np.sqrt(phase_data['acc_x']**2 + phase_data['acc_y']**2 + phase_data['acc_z']**2)\n                features[f'{phase.lower()}_acc_mag_mean'] = acc_mag.mean()\n                features[f'{phase.lower()}_acc_mag_std'] = acc_mag.std()\n                features[f'{phase.lower()}_acc_mag_max'] = acc_mag.max()\n                features[f'{phase.lower()}_acc_mag_range'] = acc_mag.max() - acc_mag.min()\n                \n                # Rotation magnitude\n                rot_mag = np.sqrt(phase_data['rot_w']**2 + phase_data['rot_x']**2 + \n                                phase_data['rot_y']**2 + phase_data['rot_z']**2)\n                features[f'{phase.lower()}_rot_mag_mean'] = rot_mag.mean()\n                features[f'{phase.lower()}_rot_mag_std'] = rot_mag.std()\n                \n                # Individual axis statistics\n                for axis in ['acc_x', 'acc_y', 'acc_z', 'rot_w', 'rot_x', 'rot_y', 'rot_z']:\n                    if axis in phase_data.columns:\n                        features[f'{phase.lower()}_{axis}_mean'] = phase_data[axis].mean()\n                        features[f'{phase.lower()}_{axis}_std'] = phase_data[axis].std()\n                        features[f'{phase.lower()}_{axis}_skew'] = phase_data[axis].skew()\n                        features[f'{phase.lower()}_{axis}_kurtosis'] = phase_data[axis].kurtosis()\n        else:\n            # Fill with zeros if phase is missing\n            features[f'{phase.lower()}_length'] = 0\n            features[f'{phase.lower()}_ratio'] = 0\n    \n    # Full sequence features\n    if all(col in seq_data.columns for col in imu_cols):\n        # Overall acceleration patterns\n        acc_mag_full = np.sqrt(seq_data['acc_x']**2 + seq_data['acc_y']**2 + seq_data['acc_z']**2)\n        features['full_acc_mag_mean'] = acc_mag_full.mean()\n        features['full_acc_mag_std'] = acc_mag_full.std()\n        features['full_acc_mag_cv'] = acc_mag_full.std() / (acc_mag_full.mean() + 1e-8)\n        \n        # Velocity estimation (numerical derivative)\n        vel_x = np.diff(seq_data['acc_x'])\n        vel_y = np.diff(seq_data['acc_y'])\n        vel_z = np.diff(seq_data['acc_z'])\n        vel_mag = np.sqrt(vel_x**2 + vel_y**2 + vel_z**2)\n        features['velocity_mean'] = vel_mag.mean()\n        features['velocity_std'] = vel_mag.std()\n        \n        # Jerk (derivative of acceleration)\n        jerk_x = np.diff(vel_x) if len(vel_x) > 1 else [0]\n        jerk_y = np.diff(vel_y) if len(vel_y) > 1 else [0]\n        jerk_z = np.diff(vel_z) if len(vel_z) > 1 else [0]\n        jerk_mag = np.sqrt(np.array(jerk_x)**2 + np.array(jerk_y)**2 + np.array(jerk_z)**2)\n        features['jerk_mean'] = jerk_mag.mean()\n        features['jerk_std'] = jerk_mag.std()\n        \n        # Zero crossings (motion pattern indicators)\n        features['acc_x_zero_crossings'] = np.sum(np.diff(np.sign(seq_data['acc_x'])) != 0)\n        features['acc_y_zero_crossings'] = np.sum(np.diff(np.sign(seq_data['acc_y'])) != 0)\n        features['acc_z_zero_crossings'] = np.sum(np.diff(np.sign(seq_data['acc_z'])) != 0)\n        \n        # Energy features\n        features['acc_energy'] = np.sum(acc_mag_full**2)\n        features['rot_energy'] = np.sum((seq_data['rot_w']**2 + seq_data['rot_x']**2 + \n                                       seq_data['rot_y']**2 + seq_data['rot_z']**2))\n        \n        # Frequency domain features (simplified FFT)\n        if len(acc_mag_full) > 4:\n            fft_vals = np.abs(fft(acc_mag_full))\n            features['fft_dominant_freq'] = np.argmax(fft_vals[1:len(fft_vals)//2]) + 1\n            features['fft_energy'] = np.sum(fft_vals**2)\n            features['fft_peak_power'] = np.max(fft_vals)\n    \n    # Thermopile features\n    if thm_cols:\n        thm_available_cols = [col for col in thm_cols if col in seq_data.columns]\n        if thm_available_cols:\n            thm_data = seq_data[thm_available_cols].dropna()\n            if not thm_data.empty:\n                features['thm_mean_temp'] = thm_data.mean().mean()\n                features['thm_std_temp'] = thm_data.std().mean()\n                features['thm_temp_range'] = thm_data.max().max() - thm_data.min().min()\n                features['thm_temp_gradient'] = (thm_data.iloc[-1] - thm_data.iloc[0]).abs().mean()\n    \n    # Time-of-flight features\n    if tof_cols:\n        tof_available_cols = [col for col in tof_cols if col in seq_data.columns]\n        if tof_available_cols:\n            tof_data = seq_data[tof_available_cols].replace(-1, np.nan)\n            tof_valid = tof_data.dropna()\n            if not tof_valid.empty:\n                features['tof_mean_distance'] = tof_valid.mean().mean()\n                features['tof_std_distance'] = tof_valid.std().mean()\n                features['tof_valid_ratio'] = tof_valid.notna().sum().sum() / (len(tof_data) * len(tof_available_cols))\n                \n                # Spatial features for ToF grids\n                for sensor_num in ['1', '2', '3', '4', '5']:\n                    sensor_cols = [col for col in tof_available_cols if f'tof_{sensor_num}_' in col]\n                    if len(sensor_cols) == 64:  # Full 8x8 grid\n                        sensor_data = tof_data[sensor_cols].replace(-1, np.nan)\n                        if not sensor_data.empty:\n                            # Reshape to 8x8 grid for spatial analysis\n                            grid_data = sensor_data.values.reshape(-1, 8, 8)\n                            valid_grids = ~np.isnan(grid_data).all(axis=(1,2))\n                            if valid_grids.any():\n                                valid_grid = grid_data[valid_grids][0]  # Take first valid grid\n                                if not np.isnan(valid_grid).all():\n                                    features[f'tof_{sensor_num}_center_mass_x'] = np.nanmean(np.arange(8)[None, :] * valid_grid) / np.nanmean(valid_grid)\n                                    features[f'tof_{sensor_num}_center_mass_y'] = np.nanmean(np.arange(8)[:, None] * valid_grid) / np.nanmean(valid_grid)\n                                    features[f'tof_{sensor_num}_spatial_variance'] = np.nanvar(valid_grid)\n    \n    return features\n\n# Extract features from sample sequences for analysis\nprint(\"Extracting features from sample sequences...\")\nsample_seq_ids = target_data['sequence_id'].unique()[:50]  # Sample 50 sequences\nfeature_data = []\n\n# Debug: Check what phases are available\nprint(\"Debug: Checking available phases in data...\")\navailable_phases = target_data['behavior'].unique()\nprint(f\"Available phases: {available_phases}\")\n\nfor seq_id in sample_seq_ids:\n    try:\n        features = extract_sequence_features(target_data, seq_id)\n        features['sequence_id'] = seq_id\n        feature_data.append(features)\n    except Exception as e:\n        print(f\"Error processing sequence {seq_id}: {e}\")\n\nfeatures_df = pd.DataFrame(feature_data)\nprint(f\"Extracted {len(features_df)} feature sets with {len(features_df.columns)} features each\")\n\n# Debug: Check what feature columns were actually created\nprint(\"\\nDebug: Feature columns created:\")\nfeature_cols = [col for col in features_df.columns if col != 'sequence_id' and col != 'gesture']\nprint(f\"Total feature columns: {len(feature_cols)}\")\nphase_feature_cols = [col for col in feature_cols if any(phase in col for phase in ['transition', 'pause', 'gesture'])]\nprint(f\"Phase-specific feature columns: {len(phase_feature_cols)}\")\nif phase_feature_cols:\n    print(\"Sample phase features:\", phase_feature_cols[:10])\nelse:\n    print(\"No phase-specific features found!\")\n\n# Analyze feature importance and distributions\nprint(\"\\nFeature analysis:\")\nnumeric_features = features_df.select_dtypes(include=[np.number]).columns\nprint(f\"Numeric features: {len(numeric_features)}\")\n\n# Plot feature distributions by gesture\ngesture_feature_analysis = features_df.groupby('gesture')[numeric_features].mean()\nprint(\"\\nTop varying features across gestures:\")\nfeature_variance = gesture_feature_analysis.var(axis=0).sort_values(ascending=False)\nprint(feature_variance.head(10))\n\n# Visualize top discriminative features\nfig, axes = plt.subplots(3, 3, figsize=(20, 15))\naxes = axes.flatten()\n\ntop_features = feature_variance.head(9).index\n\nfor i, feature in enumerate(top_features):\n    if i < len(axes):\n        # Box plot by gesture\n        gesture_data = []\n        gesture_labels = []\n        for gesture in features_df['gesture'].unique():\n            if pd.notna(features_df[features_df['gesture'] == gesture][feature]).any():\n                gesture_data.append(features_df[features_df['gesture'] == gesture][feature].dropna())\n                gesture_labels.append(gesture[:15])  # Truncate long names\n        \n        axes[i].boxplot(gesture_data, labels=gesture_labels)\n        axes[i].set_title(f'{feature}', fontsize=10)\n        axes[i].tick_params(axis='x', rotation=45, labelsize=8)\n        axes[i].grid(True, alpha=0.3)\n\nplt.tight_layout()\nplt.show()\n\n# Phase-specific analysis\nprint(\"\\nPhase-specific feature analysis:\")\n\n# First, let's check what phase features actually exist\nactual_phase_features = []\nfor phase in ['transition', 'pause', 'gesture']:\n    phase_cols = [col for col in features_df.columns if col.startswith(phase)]\n    if phase_cols:\n        actual_phase_features.extend(phase_cols)\n        print(f\"Found {len(phase_cols)} features for {phase} phase\")\n    else:\n        print(f\"No features found for {phase} phase\")\n\nif actual_phase_features:\n    # Use actual available phase features\n    available_metrics = set()\n    for col in actual_phase_features:\n        parts = col.split('_')\n        if len(parts) >= 2:\n            metric = '_'.join(parts[1:])  # Everything after phase name\n            available_metrics.add(metric)\n    \n    available_metrics = list(available_metrics)[:5]  # Take first 5 metrics\n    phase_features = ['transition', 'pause', 'gesture']\n    \n    # Create a smaller plot with actual features\n    n_metrics = min(3, len(available_metrics))\n    n_phases = 3\n    \n    if n_metrics > 0:\n        fig, axes = plt.subplots(n_metrics, n_phases, figsize=(18, 6*n_metrics))\n        if n_metrics == 1:\n            axes = axes.reshape(1, -1)\n        \n        for i, metric in enumerate(available_metrics[:n_metrics]):\n            for j, phase in enumerate(phase_features):\n                feature_name = f'{phase}_{metric}'\n                if feature_name in features_df.columns:\n                    # Distribution by gesture\n                    for k, gesture in enumerate(features_df['gesture'].unique()[:5]):  # Top 5 gestures\n                        gesture_data = features_df[features_df['gesture'] == gesture][feature_name].dropna()\n                        if not gesture_data.empty:\n                            axes[i,j].hist(gesture_data, alpha=0.6, label=gesture[:10], bins=10)\n                    \n                    axes[i,j].set_title(f'{feature_name}', fontsize=10)\n                    axes[i,j].legend(fontsize=8)\n                    axes[i,j].grid(True, alpha=0.3)\n                else:\n                    axes[i,j].text(0.5, 0.5, f'Feature {feature_name} not found', \n                                  ha='center', va='center', transform=axes[i,j].transAxes)\n                    axes[i,j].set_title(f'{feature_name} (N/A)', fontsize=10)\n        \n        plt.tight_layout()\n        plt.show()\n    else:\n        print(\"No phase-specific metrics found to plot\")\nelse:\n    print(\"No phase-specific features were generated. This might indicate an issue with:\")\n    print(\"1. Phase detection in the sequences\")\n    print(\"2. Feature extraction function\")\n    print(\"3. Data structure\")\n    \n    # Alternative: Plot available features instead\n    print(\"\\nPlotting available features instead...\")\n    available_features = [col for col in features_df.columns if col not in ['sequence_id', 'gesture', 'sequence_length']]\n    if len(available_features) >= 9:\n        fig, axes = plt.subplots(3, 3, figsize=(18, 15))\n        axes = axes.flatten()\n        \n        for i, feature in enumerate(available_features[:9]):\n            # Distribution by gesture\n            for gesture in features_df['gesture'].unique()[:5]:  # Top 5 gestures\n                gesture_data = features_df[features_df['gesture'] == gesture][feature].dropna()\n                if not gesture_data.empty:\n                    axes[i].hist(gesture_data, alpha=0.6, label=gesture[:10], bins=10)\n            \n            axes[i].set_title(f'{feature}', fontsize=10)\n            axes[i].legend(fontsize=8)\n            axes[i].grid(True, alpha=0.3)\n        \n        plt.tight_layout()\n        plt.show()\n\n# Correlation analysis of engineered features\nprint(\"\\nCorrelation analysis of engineered features:\")\ncorrelation_features = [col for col in numeric_features if 'gesture_' in col or 'full_' in col][:20]\nif correlation_features:\n    corr_matrix = features_df[correlation_features].corr()\n    \n    plt.figure(figsize=(12, 10))\n    sns.heatmap(corr_matrix, annot=True, cmap='coolwarm', center=0, fmt='.2f', \n                square=True, cbar_kws={'shrink': 0.8})\n    plt.title('Correlation Matrix of Key Engineered Features')\n    plt.xticks(rotation=45, ha='right')\n    plt.yticks(rotation=0)\n    plt.tight_layout()\n    plt.show()\n\n# Feature importance simulation\nprint(\"\\nFeature discrimination analysis:\")\nfrom scipy.stats import f_oneway\n\n# Perform ANOVA to find most discriminative features\ndiscriminative_features = []\nfor feature in numeric_features:\n    if feature in features_df.columns:\n        groups = []\n        for gesture in features_df['gesture'].unique():\n            group_data = features_df[features_df['gesture'] == gesture][feature].dropna()\n            if len(group_data) > 1:\n                groups.append(group_data)\n        \n        if len(groups) >= 2:\n            try:\n                f_stat, p_value = f_oneway(*groups)\n                if not np.isnan(f_stat):\n                    discriminative_features.append((feature, f_stat, p_value))\n            except:\n                pass\n\n# Sort by F-statistic\ndiscriminative_features.sort(key=lambda x: x[1], reverse=True)\nprint(\"\\nTop 15 most discriminative features (by ANOVA F-statistic):\")\nfor i, (feature, f_stat, p_val) in enumerate(discriminative_features[:15]):\n    print(f\"{i+1:2d}. {feature:<30} F={f_stat:.3f}, p={p_val:.6f}\")\n\n# Visualize top discriminative features\nif len(discriminative_features) >= 6:\n    fig, axes = plt.subplots(2, 3, figsize=(18, 12))\n    axes = axes.flatten()\n    \n    for i, (feature, f_stat, p_val) in enumerate(discriminative_features[:6]):\n        gesture_means = features_df.groupby('gesture')[feature].mean().sort_values(ascending=False)\n        \n        axes[i].bar(range(len(gesture_means)), gesture_means.values, alpha=0.7, edgecolor='black')\n        axes[i].set_xticks(range(len(gesture_means)))\n        axes[i].set_xticklabels([g[:10] for g in gesture_means.index], rotation=45, fontsize=8)\n        axes[i].set_title(f'{feature}\\nF={f_stat:.2f}, p={p_val:.4f}', fontsize=10)\n        axes[i].grid(True, alpha=0.3)\n    \n    plt.tight_layout()\n    plt.show()\n\nprint(\"\\n\" + \"=\"*50)\nprint(\"EDA COMPLETE\")\nprint(\"=\"*50)\nprint(\"\\nKey findings summary:\")\nprint(\"1. Dataset contains both Target (BFRB) and Non-target gestures\")\nprint(\"2. Each sequence has 3 phases: Transition, Pause, Gesture\")\nprint(\"3. Multiple sensor types: IMU (always), Thermopile, Time-of-flight\")\nprint(\"4. Missing data patterns vary by sensor type\")\nprint(\"5. Subject demographics provide additional context\")\nprint(\"6. Gesture patterns show distinct signatures in sensor data\")\nprint(\"7. Consider phase-specific feature extraction for better performance\")\nprint(\"8. Engineered features show strong discriminative power between gestures\")\nprint(\"9. Gesture-phase features appear most informative for classification\")\nprint(\"10. Acceleration and rotation magnitude features are key discriminators\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T18:03:14.513328Z","iopub.execute_input":"2025-05-30T18:03:14.513899Z","iopub.status.idle":"2025-05-30T18:04:12.207095Z","shell.execute_reply.started":"2025-05-30T18:03:14.513866Z","shell.execute_reply":"2025-05-30T18:04:12.205159Z"},"_kg_hide-input":true},"outputs":[],"execution_count":null}]}