OSCR

EEG-Based Supported Diagnosis of ADHD Using Subject-Specific HMMs and Stationary RKHS Embeddings.

Code ↔ Paper

12 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 12 matches · 2 of them tie a paragraph to a whole file without given lines: weak matches, whose lines are not tinted
  1. [1] § 2. Materials and Methods › 2.3. Preprocessing ↔ baseline/baseline.py, lines 131 to 143 · score 0.83 · 13–30 Hz, 0.5–4 Hz, 8–13 Hz, 4–8 Hz, theta, delta · Report a mismatch
  2. [2] § 3. Results › 3.1. Performance Analysis on the Synthetic Dataset ↔ single held out test/rkhs-knn-svm-eeg-dataset.ipynb, lines 1037 to 1185 · score 0.77 · full training, RKHS kernel, AUC ROC, Optuna, cross validation, Matthews · Report a mismatch
  3. [3] § 3. Results › 3.3. Performance Analysis on the Real EEG Dataset ↔ baseline/baseline.py, lines 1 to 61 · score 0.73 · Alim, Imtiaz, reimplemented, Hjorth, baseline, theta · Report a mismatch
  4. [4] § 2. Materials and Methods › 2.4. Hidden Markov Model-Based Time-Series Representation ↔ training-hmm.ipynb, lines 649 to 745 · score 0.71 · Baum Welch, log likelihood, EM, iterations, fit, training · Report a mismatch
  5. [5] § 2. Materials and Methods › 2.1. Synthetic EEG Benchmark Signals ↔ Synthetic_EEG_data_generation.ipynb, the whole file · a weak match · score 0.66 · frequency bands, EEG signals, healthy, noise, simulated, synthetic · Report a mismatch
  6. [6] § 3. Results › 3.2. Hyperparameter Tuning ↔ single held out test/rkhs-knn-svm-eeg-dataset.ipynb, lines 1037 to 1185 · score 0.65 · kernel transformation, AUC ROC, confusion matrices, Matthews, balanced accuracy, bandwidth · Report a mismatch
  7. [7] § 2. Materials and Methods ↔ single held out test/rkhs-knn-svm-eeg-dataset.ipynb, lines 832 to 908 · score 0.62 · HMM GMM model, distance computation, precomputed kernel, RKHS, EEG · Report a mismatch
  8. [8] § 2. Materials and Methods › 2.1. Synthetic EEG Benchmark Signals ↔ baseline/baseline.py, lines 1 to 61 · score 0.61 · wave, Hz, brain, theta, phase, delta · Report a mismatch
  9. [9] § 2. Materials and Methods › 2.6. Classification with Precomputed Kernels ↔ single held out test/rkhs-knn-svm-eeg-dataset.ipynb, lines 798 to 830 · score 0.60 · positive semi definite, kernel matrix, SVM · Report a mismatch
  10. [10] § 2. Materials and Methods › 2.1. Synthetic EEG Benchmark Signals ↔ Synthetic_EEG_data_generation.ipynb, the whole file · a weak match · score 0.57 · wave, uniform, noise, theta, delta, simulated · Report a mismatch
  11. [11] § 3. Results › 3.2. Hyperparameter Tuning ↔ training-hmm.ipynb, lines 1143 to 1229 · score 0.52 · Log likelihood, HMM convergence, iteration, errors, fit · Report a mismatch
  12. [12] § 2. Materials and Methods › 2.6. Classification with Precomputed Kernels ↔ nested CV/PPK_with_knn_and_svm.ipynb, lines 293 to 317 · score 0.51 · Gram matrix, unnecessary, majority, formulation, SVM, distances · Report a mismatch

Paper

Loaded from Europe PMC by your browser; OSCR does not store it: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Jupyter notebook · 1,221 lines · 47 KB · MIT · 4 matches

  1. # %%
  2. # --- Install required packages ---
  3. !pip install -q tensorflow tensorflow-probability scikit-learn matplotlib pandas optuna
  5. # %% [markdown]
  6. # # RKHS with precomputed distance and KNN Classification
  8. # %%
  9. import numpy as np
  10. import pandas as pd
  11. from sklearn.neighbors import KNeighborsClassifier
  12. from sklearn.model_selection import StratifiedKFold, train_test_split
  13. from sklearn.metrics import (
  14. classification_report, accuracy_score, confusion_matrix,
  15. balanced_accuracy_score, roc_auc_score, matthews_corrcoef,
  16. roc_curve, precision_recall_curve, auc
  17. )
  18. import matplotlib.pyplot as plt
  19. import seaborn as sns
  20. import optuna
  21. from optuna.samplers import TPESampler
  22. import warnings
  23. import scipy.linalg
  24. warnings.filterwarnings('ignore')
  26. # -----------------------------
  27. # Stationary distribution utils
  28. # -----------------------------
  29. def compute_stationary_distribution(A, max_iter=1000, rel_tol=1e-6, abs_tol=1e-10):
  30. A = np.asarray(A, dtype=float)
  31. n = A.shape[0]
  33. # Row-stochastic guard (and all-zero row fix)
  34. rs = A.sum(axis=1, keepdims=True)
  35. bad = (rs[:, 0] == 0)
  36. if bad.any():
  37. A[bad] = 1.0 / n
  38. rs = A.sum(axis=1, keepdims=True)
  39. A = A / rs
  41. pi = np.ones(n) / n
  42. for _ in range(max_iter):
  43. pi_new = pi @ A
  44. if np.linalg.norm(pi_new - pi, 1) <= rel_tol * np.linalg.norm(pi, 1) + abs_tol:
  45. pi = pi_new
  46. break
  47. pi = pi_new
  49. pi = np.maximum(pi, 0.0)
  50. s = pi.sum()
  51. return pi / (s if s > 0 else 1.0)
  53. # -----------------------------------------
  54. # Exact multivariate Gaussian-Gaussian RBF
  55. # -----------------------------------------
  56. def k_gauss_of_gaussians(mu1, S1, mu2, S2, ell, eps=1e-9):
  57. mu1 = np.asarray(mu1); mu2 = np.asarray(mu2)
  58. S1 = np.asarray(S1); S2 = np.asarray(S2)
  59. D = mu1.shape[0]
  61. S1 = S1 + eps * np.eye(D)
  62. S2 = S2 + eps * np.eye(D)
  64. L_det = np.linalg.cholesky(np.eye(D) + (S1 + S2) / max(ell, eps))
  65. log_det_term = -np.sum(np.log(np.diag(L_det)))
  67. M = S1 + S2 + ell * np.eye(D)
  68. Lm = np.linalg.cholesky(M)
  69. d = mu1 - mu2
  70. y = scipy.linalg.solve_triangular(Lm, d, lower=True)
  71. quad = -0.5 * float(y @ y)
  73. return float(np.exp(log_det_term + quad))
  75. # -------------------------------------------------------
  76. # CSV extraction, with auto-detected feature dimensionality
  77. # -------------------------------------------------------
  78. def infer_n_features(df, n_states, n_components):
  79. f_idxs = []
  80. prefix = f"gmm_mean_0_0_f"
  81. for c in df.columns:
  82. if c.startswith(prefix):
  83. try:
  84. f_idxs.append(int(c.split('f')[-1]))
  85. except Exception:
  86. pass
  87. return (max(f_idxs) + 1) if f_idxs else 1
  89. def extract_hmm_parameters_from_csv(df, n_states=3, n_components=3, n_features=None):
  90. if n_features is None:
  91. n_features = infer_n_features(df, n_states, n_components)
  93. models = []
  94. print(f"Extracting HMM parameters from {len(df)} models...")
  95. print(f"Config inferred: n_states={n_states}, n_components={n_components}, n_features={n_features}")
  97. has_pi_s = all([f"stationary_pi_{i}" in df.columns for i in range(n_states)])
  98. has_w_eff = f"w_eff_0_0" in df.columns
  100. for idx, row in df.iterrows():
  101. try:
  102. # Transition matrix A
  103. A = np.zeros((n_states, n_states), dtype=float)
  104. for i in range(n_states):
  105. for j in range(n_states):
  106. A[i, j] = row[f"A_{i}{j}"]
  107. rs = A.sum(axis=1, keepdims=True); rs[rs == 0] = 1
  108. A = A / rs
  110. # stationary pi
  111. if has_pi_s:
  112. pi_s = np.array([row[f"stationary_pi_{i}"] for i in range(n_states)], dtype=float)
  113. if (not np.isfinite(pi_s).all()) or (pi_s.sum() <= 1e-12) or (pi_s < -1e-12).any():
  114. pi_s = compute_stationary_distribution(A)
  115. else:
  116. pi_s = np.maximum(pi_s, 0.0)
  117. s = pi_s.sum()
  118. pi_s = pi_s / (s if s > 0 else 1.0)
  119. else:
  120. pi_s = compute_stationary_distribution(A)
  122. # mixture weights alpha and means/covs
  123. alpha = np.zeros((n_states, n_components), dtype=float)
  124. mu = np.zeros((n_states, n_components, n_features), dtype=float)
  125. sigma = np.zeros((n_states, n_components, n_features, n_features), dtype=float)
  127. for i in range(n_states):
  128. for k in range(n_components):
  129. if has_w_eff:
  130. alpha[i, k] = row.get(f"gmm_weight_{i}_{k}", 0.0)
  131. else:
  132. alpha[i, k] = row[f"gmm_weight_{i}_{k}"]
  134. for f in range(n_features):
  135. mu[i, k, f] = row[f"gmm_mean_{i}_{k}_f{f}"]
  137. Si = np.zeros((n_features, n_features), dtype=float)
  138. for f1 in range(n_features):
  139. for f2 in range(f1, n_features):
  140. Si[f1, f2] = row[f"gmm_cov_{i}_{k}_f{f1}f{f2}"]
  141. if f1 != f2:
  142. Si[f2, f1] = Si[f1, f2]
  143. Si += 1e-6 * np.eye(n_features)
  144. try:
  145. np.linalg.cholesky(Si)
  146. except np.linalg.LinAlgError:
  147. Si = np.eye(n_features) * 0.1
  148. sigma[i, k] = Si
  150. if not has_w_eff:
  151. alpha = np.maximum(alpha, 0.0)
  152. sums = alpha.sum(axis=1, keepdims=True); sums[sums == 0] = 1.0
  153. alpha = alpha / sums
  155. model = {
  156. 'pi_stationary': pi_s,
  157. 'alpha': alpha,
  158. 'mu': mu,
  159. 'sigma': sigma,
  160. 'transition_matrix': A,
  161. 'n_states': n_states,
  162. 'n_components': n_components,
  163. 'n_features': n_features,
  164. 'model_id': idx,
  165. 'dataset_type': row.get('dataset_type', 'unknown'),
  166. 'series_id': row.get('series_id', idx)
  167. }
  169. if has_w_eff:
  170. w_eff = np.array(
  171. [row[f"w_eff_{i}_{k}"] for i in range(n_states) for k in range(n_components)],
  172. dtype=float
  173. )
  174. model['w_eff'] = w_eff
  176. models.append(model)
  178. except Exception as e:
  179. print(f"Error processing model {idx}: {str(e)}")
  180. continue
  182. print(f"Successfully extracted parameters for {len(models)} models")
  183. return models
  185. # --------------------------
  186. # Collapse to component lists
  187. # --------------------------
  188. def model_to_components(model):
  189. K = model['n_states']; M = model['n_components']; D = model['n_features']
  190. pi_s = model['pi_stationary']
  191. alpha = model['alpha']
  192. mu = model['mu']
  193. sigma = model['sigma']
  195. if 'w_eff' in model:
  196. w = model['w_eff'].copy()
  197. else:
  198. w = np.array([pi_s[i] * alpha[i, k] for i in range(K) for k in range(M)], dtype=float)
  200. MU = np.array([mu[i, k, :] for i in range(K) for k in range(M)], dtype=float)
  201. SIG = np.array([sigma[i, k, :, :] for i in range(K) for k in range(M)], dtype=float)
  202. return w, MU, SIG
  204. # -------------------------------------
  205. # Distance (MMD²) and matrix computation
  206. # -------------------------------------
  207. def enhanced_rkhs_distance_multivariate(modelP, modelQ, sigma_kernel):
  208. wP, MUP, SIGP = model_to_components(modelP)
  209. wQ, MUQ, SIGQ = model_to_components(modelQ)
  211. term_PP = 0.0
  212. for i in range(len(wP)):
  213. for j in range(len(wP)):
  214. term_PP += wP[i] * wP[j] * k_gauss_of_gaussians(MUP[i], SIGP[i], MUP[j], SIGP[j], sigma_kernel)
  216. term_QQ = 0.0
  217. for i in range(len(wQ)):
  218. for j in range(len(wQ)):
  219. term_QQ += wQ[i] * wQ[j] * k_gauss_of_gaussians(MUQ[i], SIGQ[i], MUQ[j], SIGQ[j], sigma_kernel)
  221. term_PQ = 0.0
  222. for i in range(len(wP)):
  223. for j in range(len(wQ)):
  224. term_PQ += wP[i] * wQ[j] * k_gauss_of_gaussians(MUP[i], SIGP[i], MUQ[j], SIGQ[j], sigma_kernel)
  226. dist2 = term_PP + term_QQ - 2.0 * term_PQ
  227. return max(float(dist2), 0.0)
  229. def compute_distance_matrix(models, sigma_kernel):
  230. n = len(models)
  231. D = np.zeros((n, n), dtype=float)
  233. print(f"Computing {n}x{n} RKHS distance matrix (σ={sigma_kernel:.6g})...")
  234. comps = [model_to_components(m) for m in models]
  236. def dist2_from_comps(cP, cQ):
  237. wP, MUP, SIGP = cP
  238. wQ, MUQ, SIGQ = cQ
  239. term_PP = 0.0
  240. for i in range(len(wP)):
  241. for j in range(len(wP)):
  242. term_PP += wP[i] * wP[j] * k_gauss_of_gaussians(MUP[i], SIGP[i], MUP[j], SIGP[j], sigma_kernel)
  243. term_QQ = 0.0
  244. for i in range(len(wQ)):
  245. for j in range(len(wQ)):
  246. term_QQ += wQ[i] * wQ[j] * k_gauss_of_gaussians(MUQ[i], SIGQ[i], MUQ[j], SIGQ[j], sigma_kernel)
  247. term_PQ = 0.0
  248. for i in range(len(wP)):
  249. for j in range(len(wQ)):
  250. term_PQ += wP[i] * wQ[j] * k_gauss_of_gaussians(MUP[i], SIGP[i], MUQ[j], SIGQ[j], sigma_kernel)
  251. return max(float(term_PP + term_QQ - 2.0 * term_PQ), 0.0)
  253. for i in range(n):
  254. D[i, i] = 0.0
  255. for j in range(i + 1, n):
  256. d2 = dist2_from_comps(comps[i], comps[j])
  257. D[i, j] = D[j, i] = d2
  258. if (i + 1) % 20 == 0:
  259. print(f" Processed {i + 1}/{n} rows")
  261. print(f"Done. Stats: min={D.min():.6f}, max={D.max():.6f}, mean={D.mean():.6f}")
  262. return D
  264. # ----------------------------
  265. # Cross-validation (precompute)
  266. # ----------------------------
  267. def kfold_cross_validation(models, labels, sigma_kernel, k_neighbors, n_splits=5, random_state=42):
  268. skf = StratifiedKFold(n_splits=n_splits, shuffle=True, random_state=random_state)
  269. full_D = compute_distance_matrix(models, sigma_kernel)
  271. fold_results = {
  272. 'accuracy': [], 'balanced_accuracy': [], 'auc_roc': [], 'mcc': [],
  273. 'predictions': [], 'probabilities': [], 'true_labels': [],
  274. 'confusion_matrices': []
  275. }
  277. print(f"Performing {n_splits}-fold cross-validation...")
  278. for fold, (train_idx, val_idx) in enumerate(skf.split(models, labels)):
  279. print(f" Fold {fold + 1}/{n_splits}")
  281. train_D = full_D[np.ix_(train_idx, train_idx)]
  282. val_train_D = full_D[np.ix_(val_idx, train_idx)]
  284. knn = KNeighborsClassifier(n_neighbors=k_neighbors, metric='precomputed', weights='distance')
  285. knn.fit(train_D, labels[train_idx])
  287. val_pred = knn.predict(val_train_D)
  288. val_prob = knn.predict_proba(val_train_D)[:, 1]
  290. fold_results['predictions'].extend(val_pred.tolist())
  291. fold_results['probabilities'].extend(val_prob.tolist())
  292. fold_results['true_labels'].extend(labels[val_idx].tolist())
  294. acc = accuracy_score(labels[val_idx], val_pred)
  295. bacc = balanced_accuracy_score(labels[val_idx], val_pred)
  296. auc_ = roc_auc_score(labels[val_idx], val_prob)
  297. mcc = matthews_corrcoef(labels[val_idx], val_pred)
  298. cm = confusion_matrix(labels[val_idx], val_pred)
  300. fold_results['accuracy'].append(acc)
  301. fold_results['balanced_accuracy'].append(bacc)
  302. fold_results['auc_roc'].append(auc_)
  303. fold_results['mcc'].append(mcc)
  304. fold_results['confusion_matrices'].append(cm)
  306. print(f" Acc={acc:.4f}, BalAcc={bacc:.4f}, AUC={auc_:.4f}, MCC={mcc:.4f}")
  308. for k in ['accuracy', 'balanced_accuracy', 'auc_roc', 'mcc']:
  309. fold_results[k] = np.array(fold_results[k])
  310. fold_results['predictions'] = np.array(fold_results['predictions'])
  311. fold_results['probabilities'] = np.array(fold_results['probabilities'])
  312. fold_results['true_labels'] = np.array(fold_results['true_labels'])
  313. return fold_results
  315. # -----------------------------
  316. # Optuna (precompute per sigma too)
  317. # -----------------------------
  318. def hyperparameter_optimization(models, labels, sigma_kernel_range=(1e-6, 1.0),
  319. k_neighbors_range=(1, 15), n_trials=50,
  320. n_splits=5, timeout=1800):
  321. print("Starting hyperparameter optimization...")
  323. min_k, max_k = k_neighbors_range
  324. max_k = min(max_k, len(models) - 1)
  325. valid_k_values = [k for k in range(min_k, max_k + 1) if k % 2 == 1]
  327. print(f"Search space:")
  328. print(f"- sigma_kernel: log-uniform from {sigma_kernel_range[0]:.1e} to {sigma_kernel_range[1]:.1e}")
  329. print(f"- k_neighbors: {valid_k_values}")
  330. print(f"- n_trials: {n_trials}, CV folds: {n_splits}")
  332. def objective(trial):
  333. sigma_kernel = trial.suggest_float('sigma_kernel', sigma_kernel_range[0], sigma_kernel_range[1], log=True)
  334. k_neighbors = trial.suggest_categorical('k_neighbors', valid_k_values)
  336. try:
  337. full_D = compute_distance_matrix(models, sigma_kernel)
  338. skf = StratifiedKFold(n_splits=n_splits, shuffle=True, random_state=42)
  339. scores = []
  340. for train_idx, val_idx in skf.split(models, labels):
  341. train_D = full_D[np.ix_(train_idx, train_idx)]
  342. val_train_D = full_D[np.ix_(val_idx, train_idx)]
  343. knn = KNeighborsClassifier(n_neighbors=k_neighbors, metric='precomputed', weights='distance')
  344. knn.fit(train_D, labels[train_idx])
  345. pred = knn.predict(val_train_D)
  346. scores.append(balanced_accuracy_score(labels[val_idx], pred))
  347. return float(np.mean(scores))
  348. except Exception as e:
  349. print(f"Trial failed: {e}")
  350. return 0.0
  352. study = optuna.create_study(direction='maximize', sampler=TPESampler(seed=42),
  353. study_name="RKHS_KNN_Optimization")
  354. study.optimize(objective, n_trials=n_trials, timeout=timeout)
  356. print("\nOptimization completed!")
  357. print(f"Best balanced accuracy: {study.best_value:.4f}")
  358. print(f"Best parameters: {study.best_params}")
  359. return study
  361. # -----------------
  362. # Evaluation plots
  363. # -----------------
  364. def create_evaluation_plots(fold_results, class_names=['Control', 'ADHD']):
  365. fig, axes = plt.subplots(2, 3, figsize=(18, 12))
  367. overall_cm = confusion_matrix(fold_results['true_labels'], fold_results['predictions'])
  368. sns.heatmap(overall_cm, annot=True, fmt='d', cmap='Blues',
  369. xticklabels=class_names, yticklabels=class_names, ax=axes[0,0])
  370. axes[0,0].set_title('Overall Confusion Matrix'); axes[0,0].set_xlabel('Predicted'); axes[0,0].set_ylabel('Actual')
  372. fpr, tpr, _ = roc_curve(fold_results['true_labels'], fold_results['probabilities'])
  373. auc_score = auc(fpr, tpr)
  374. axes[0,1].plot(fpr, tpr, lw=2, label=f'ROC (AUC={auc_score:.3f})')
  375. axes[0,1].plot([0, 1], [0, 1], lw=2, linestyle='--', label='Random')
  376. axes[0,1].set_xlabel('FPR'); axes[0,1].set_ylabel('TPR'); axes[0,1].set_title('ROC Curve'); axes[0,1].legend(); axes[0,1].grid(True, alpha=0.3)
  378. precision, recall, _ = precision_recall_curve(fold_results['true_labels'], fold_results['probabilities'])
  379. pr_auc = auc(recall, precision)
  380. axes[0,2].plot(recall, precision, lw=2, label=f'PR (AUC={pr_auc:.3f})')
  381. axes[0,2].set_xlabel('Recall'); axes[0,2].set_ylabel('Precision'); axes[0,2].set_title('Precision-Recall'); axes[0,2].legend(); axes[0,2].grid(True, alpha=0.3)
  383. metrics_data = [fold_results['accuracy'], fold_results['balanced_accuracy'], fold_results['auc_roc'], fold_results['mcc']]
  384. axes[1,0].boxplot(metrics_data, labels=['Acc', 'Bal.Acc', 'AUC', 'MCC'])
  385. axes[1,0].set_title('CV Metric Distributions'); axes[1,0].set_ylabel('Score'); axes[1,0].grid(True, alpha=0.3)
  387. folds = np.arange(1, len(fold_results['accuracy']) + 1)
  388. axes[1,1].plot(folds, fold_results['balanced_accuracy'], 'o-', label='Balanced Acc')
  389. axes[1,1].plot(folds, fold_results['auc_roc'], 'o-', label='AUC-ROC')
  390. axes[1,1].set_xlabel('Fold'); axes[1,1].set_ylabel('Score'); axes[1,1].set_title('Performance by Fold'); axes[1,1].legend(); axes[1,1].grid(True, alpha=0.3)
  392. class_counts = np.bincount(fold_results['true_labels'])
  393. axes[1,2].pie(class_counts, labels=class_names, autopct='%1.1f%%', startangle=90)
  394. axes[1,2].set_title('Class Distribution')
  396. plt.tight_layout(); plt.show()
  398. # -----------------------
  399. # Full analysis pipeline
  400. # -----------------------
  401. def run_complete_analysis(csv_path, test_size=0.2, random_state=42,
  402. sigma_kernel_range=(1e-6, 1.0), k_neighbors_range=(1, 15),
  403. n_trials=50, n_splits=5, timeout=1800,
  404. n_states=3, n_components=3, n_features=None):
  405. print("=== RKHS-KNN ANALYSIS FOR HMM-GMM MODELS ===\n")
  407. print("Loading CSV...")
  408. df = pd.read_csv(csv_path)
  409. df['class_label'] = (df['dataset_type'] == 'ADHD').astype(int)
  411. print(f"Dataset summary:")
  412. print(f" Total: {len(df)} | ADHD: {df['class_label'].sum()} | Control: {(1-df['class_label']).sum()}")
  413. print(f" Class balance: {df['class_label'].mean():.3f}")
  415. train_df, test_df = train_test_split(
  416. df, test_size=test_size, random_state=random_state, stratify=df['class_label']
  417. )
  419. if n_features is None:
  420. n_features = infer_n_features(train_df, n_states, n_components)
  422. print("\nExtracting HMM-GMM parameters...")
  423. train_models = extract_hmm_parameters_from_csv(train_df, n_states, n_components, n_features)
  424. train_labels = train_df['class_label'].values
  425. test_models = extract_hmm_parameters_from_csv(test_df, n_states, n_components, n_features)
  426. test_labels = test_df['class_label'].values
  428. study = hyperparameter_optimization(
  429. train_models, train_labels,
  430. sigma_kernel_range=sigma_kernel_range,
  431. k_neighbors_range=k_neighbors_range,
  432. n_trials=n_trials, n_splits=n_splits, timeout=timeout
  433. )
  434. best_params = study.best_params
  436. print("\nCross-validation with best params...")
  437. fold_results = kfold_cross_validation(
  438. train_models, train_labels,
  439. best_params['sigma_kernel'], best_params['k_neighbors'],
  440. n_splits=n_splits, random_state=random_state
  441. )
  443. print("\n" + "="*60)
  444. print("CROSS-VALIDATION RESULTS")
  445. print("="*60)
  446. print(f" σ_kernel: {best_params['sigma_kernel']:.6g}")
  447. print(f" k_neighbors: {best_params['k_neighbors']}")
  448. for name, key in [('Accuracy', 'accuracy'), ('Balanced Accuracy', 'balanced_accuracy'),
  449. ('AUC-ROC','auc_roc'), ('MCC','mcc')]:
  450. vals = fold_results[key]
  451. print(f" {name}: {vals.mean():.4f} ± {vals.std():.4f}")
  453. print("\n" + "="*60)
  454. print("FINAL TEST SET EVALUATION")
  455. print("="*60)
  456. train_D = compute_distance_matrix(train_models, best_params['sigma_kernel'])
  457. all_models = train_models + test_models
  458. full_D = compute_distance_matrix(all_models, best_params['sigma_kernel'])
  459. n_train = len(train_models)
  460. test_train_D = full_D[n_train:, :n_train]
  462. final_knn = KNeighborsClassifier(n_neighbors=best_params['k_neighbors'], metric='precomputed')
  463. final_knn.fit(train_D, train_labels)
  465. test_pred = final_knn.predict(test_train_D)
  466. test_prob = final_knn.predict_proba(test_train_D)[:, 1]
  468. test_accuracy = accuracy_score(test_labels, test_pred)
  469. test_bal_acc = balanced_accuracy_score(test_labels, test_pred)
  470. test_auc = roc_auc_score(test_labels, test_prob)
  471. test_mcc = matthews_corrcoef(test_labels, test_pred)
  472. test_cm = confusion_matrix(test_labels, test_pred)
  474. print(f" Accuracy: {test_accuracy:.4f}")
  475. print(f" Balanced Accuracy: {test_bal_acc:.4f}")
  476. print(f" AUC-ROC: {test_auc:.4f}")
  477. print(f" Matthews Corr: {test_mcc:.4f}")
  478. print("\nTest Confusion Matrix:")
  479. print(f" Control ADHD")
  480. print(f"Control {test_cm[0,0]:<7} {test_cm[0,1]:<7}")
  481. print(f"ADHD {test_cm[1,0]:<7} {test_cm[1,1]:<7}")
  483. create_evaluation_plots(fold_results)
  485. return {
  486. 'study': study,
  487. 'best_params': best_params,
  488. 'fold_results': fold_results,
  489. 'test_results': {
  490. 'accuracy': test_accuracy,
  491. 'balanced_accuracy': test_bal_acc,
  492. 'auc_roc': test_auc,
  493. 'mcc': test_mcc,
  494. 'confusion_matrix': test_cm
  495. }
  496. }
  498. if __name__ == "__main__":
  499. results = run_complete_analysis(
  500. # csv_path='.../synthetic_hmms_results.csv',
  501. csv_path='.../NxGx/hmm_results.csv',
  502. test_size=0.2,
  503. random_state=42,
  504. sigma_kernel_range=(0.001, 0.9),
  505. k_neighbors_range=(1, 11),
  506. n_trials=200,
  507. n_splits=7, # k-fold
  508. timeout=(60 * 60 * 3), # 60[sec]*60[min]*X[hour]
  509. n_states=3, #4, 5
  510. n_components=3, #4 ,5
  511. n_features=None
  512. )
  515. print("\n" + "="*60)
  516. print("ANALYSIS COMPLETE")
  517. print("="*60)
  518. print(f"Best σ_kernel: {results['best_params']['sigma_kernel']:.6g}")
  519. print(f"Best k_neighbors: {results['best_params']['k_neighbors']}")
  520. print(f"Final test balanced accuracy: {results['test_results']['balanced_accuracy']:.4f}")
  522. # %% [markdown]
  523. # # RKHS with precomputed distance and SVM Classification
  525. # %%
  526. import numpy as np
  527. import pandas as pd
  528. from sklearn.svm import SVC
  529. from sklearn.model_selection import StratifiedKFold, train_test_split
  530. from sklearn.metrics import (
  531. classification_report, accuracy_score, confusion_matrix,
  532. balanced_accuracy_score, roc_auc_score, matthews_corrcoef,
  533. roc_curve, precision_recall_curve, auc
  534. )
  535. import matplotlib.pyplot as plt
  536. import seaborn as sns
  537. import optuna
  538. from optuna.samplers import TPESampler
  539. import warnings
  540. import scipy.linalg
  541. warnings.filterwarnings('ignore')
  543. # -----------------------------
  544. # Stationary distribution utils
  545. # -----------------------------
  546. def compute_stationary_distribution(A, max_iter=1000, rel_tol=1e-6, abs_tol=1e-10):
  547. A = np.asarray(A, dtype=float)
  548. n = A.shape[0]
  550. rs = A.sum(axis=1, keepdims=True)
  551. bad = (rs[:, 0] == 0)
  552. if bad.any():
  553. A[bad] = 1.0 / n
  554. rs = A.sum(axis=1, keepdims=True)
  555. A = A / rs
  557. pi = np.ones(n) / n
  558. for _ in range(max_iter):
  559. pi_new = pi @ A
  560. if np.linalg.norm(pi_new - pi, 1) <= rel_tol * np.linalg.norm(pi, 1) + abs_tol:
  561. pi = pi_new
  562. break
  563. pi = pi_new
  565. pi = np.maximum(pi, 0.0)
  566. s = pi.sum()
  567. return pi / (s if s > 0 else 1.0)
  569. # -----------------------------------------
  570. # Exact multivariate Gaussian-Gaussian RBF
  571. # -----------------------------------------
  572. def k_gauss_of_gaussians(mu1, S1, mu2, S2, ell, eps=1e-9):
  573. mu1 = np.asarray(mu1); mu2 = np.asarray(mu2)
  574. S1 = np.asarray(S1); S2 = np.asarray(S2)
  575. D = mu1.shape[0]
  577. S1 = S1 + eps * np.eye(D)
  578. S2 = S2 + eps * np.eye(D)
  580. L_det = np.linalg.cholesky(np.eye(D) + (S1 + S2) / max(ell, eps))
  581. log_det_term = -np.sum(np.log(np.diag(L_det)))
  583. M = S1 + S2 + ell * np.eye(D)
  584. Lm = np.linalg.cholesky(M)
  585. d = mu1 - mu2
  586. y = scipy.linalg.solve_triangular(Lm, d, lower=True)
  587. quad = -0.5 * float(y @ y)
  589. return float(np.exp(log_det_term + quad))
  591. # -------------------------------------------------------
  592. # CSV extraction, with auto-detected feature dimensionality
  593. # -------------------------------------------------------
  594. def infer_n_features(df, n_states, n_components):
  595. f_idxs = []
  596. prefix = f"gmm_mean_0_0_f"
  597. for c in df.columns:
  598. if c.startswith(prefix):
  599. try:
  600. f_idxs.append(int(c.split('f')[-1]))
  601. except Exception:
  602. pass
  603. return (max(f_idxs) + 1) if f_idxs else 1
  605. def extract_hmm_parameters_from_csv(df, n_states=3, n_components=3, n_features=None):
  606. if n_features is None:
  607. n_features = infer_n_features(df, n_states, n_components)
  609. models = []
  610. print(f"Extracting HMM parameters from {len(df)} models...")
  611. print(f"Config inferred: n_states={n_states}, n_components={n_components}, n_features={n_features}")
  613. has_pi_s = all([f"stationary_pi_{i}" in df.columns for i in range(n_states)])
  614. has_w_eff = f"w_eff_0_0" in df.columns
  616. for idx, row in df.iterrows():
  617. try:
  618. A = np.zeros((n_states, n_states), dtype=float)
  619. for i in range(n_states):
  620. for j in range(n_states):
  621. A[i, j] = row[f"A_{i}{j}"]
  622. rs = A.sum(axis=1, keepdims=True); rs[rs == 0] = 1
  623. A = A / rs
  625. if has_pi_s:
  626. pi_s = np.array([row[f"stationary_pi_{i}"] for i in range(n_states)], dtype=float)
  627. if (not np.isfinite(pi_s).all()) or (pi_s.sum() <= 1e-12) or (pi_s < -1e-12).any():
  628. pi_s = compute_stationary_distribution(A)
  629. else:
  630. pi_s = np.maximum(pi_s, 0.0)
  631. s = pi_s.sum()
  632. pi_s = pi_s / (s if s > 0 else 1.0)
  633. else:
  634. pi_s = compute_stationary_distribution(A)
  636. alpha = np.zeros((n_states, n_components), dtype=float)
  637. mu = np.zeros((n_states, n_components, n_features), dtype=float)
  638. sigma = np.zeros((n_states, n_components, n_features, n_features), dtype=float)
  640. for i in range(n_states):
  641. for k in range(n_components):
  642. if has_w_eff:
  643. alpha[i, k] = row.get(f"gmm_weight_{i}_{k}", 0.0)
  644. else:
  645. alpha[i, k] = row[f"gmm_weight_{i}_{k}"]
  647. for f in range(n_features):
  648. mu[i, k, f] = row[f"gmm_mean_{i}_{k}_f{f}"]
  650. Si = np.zeros((n_features, n_features), dtype=float)
  651. for f1 in range(n_features):
  652. for f2 in range(f1, n_features):
  653. Si[f1, f2] = row[f"gmm_cov_{i}_{k}_f{f1}f{f2}"]
  654. if f1 != f2:
  655. Si[f2, f1] = Si[f1, f2]
  656. Si += 1e-6 * np.eye(n_features)
  657. try:
  658. np.linalg.cholesky(Si)
  659. except np.linalg.LinAlgError:
  660. Si = np.eye(n_features) * 0.1
  661. sigma[i, k] = Si
  663. if not has_w_eff:
  664. alpha = np.maximum(alpha, 0.0)
  665. sums = alpha.sum(axis=1, keepdims=True); sums[sums == 0] = 1.0
  666. alpha = alpha / sums
  668. model = {
  669. 'pi_stationary': pi_s,
  670. 'alpha': alpha,
  671. 'mu': mu,
  672. 'sigma': sigma,
  673. 'transition_matrix': A,
  674. 'n_states': n_states,
  675. 'n_components': n_components,
  676. 'n_features': n_features,
  677. 'model_id': idx,
  678. 'dataset_type': row.get('dataset_type', 'unknown'),
  679. 'series_id': row.get('series_id', idx)
  680. }
  682. if has_w_eff:
  683. w_eff = np.array(
  684. [row[f"w_eff_{i}_{k}"] for i in range(n_states) for k in range(n_components)],
  685. dtype=float
  686. )
  687. model['w_eff'] = w_eff
  689. models.append(model)
  691. except Exception as e:
  692. print(f"Error processing model {idx}: {str(e)}")
  693. continue
  695. print(f"Successfully extracted parameters for {len(models)} models")
  696. return models
  698. # --------------------------
  699. # Collapse to component lists
  700. # --------------------------
  701. def model_to_components(model):
  702. K = model['n_states']; M = model['n_components']; D = model['n_features']
  703. pi_s = model['pi_stationary']
  704. alpha = model['alpha']
  705. mu = model['mu']
  706. sigma = model['sigma']
  708. if 'w_eff' in model:
  709. w = model['w_eff'].copy()
  710. else:
  711. w = np.array([pi_s[i] * alpha[i, k] for i in range(K) for k in range(M)], dtype=float)
  713. MU = np.array([mu[i, k, :] for i in range(K) for k in range(M)], dtype=float)
  714. SIG = np.array([sigma[i, k, :, :] for i in range(K) for k in range(M)], dtype=float)
  715. return w, MU, SIG
  717. # -------------------------------------
  718. # Distance (MMD²) and matrix computation
  719. # -------------------------------------
  720. def enhanced_rkhs_distance_multivariate(modelP, modelQ, sigma_kernel):
  721. wP, MUP, SIGP = model_to_components(modelP)
  722. wQ, MUQ, SIGQ = model_to_components(modelQ)
  724. term_PP = 0.0
  725. for i in range(len(wP)):
  726. for j in range(len(wP)):
  727. term_PP += wP[i] * wP[j] * k_gauss_of_gaussians(MUP[i], SIGP[i], MUP[j], SIGP[j], sigma_kernel)
  729. term_QQ = 0.0
  730. for i in range(len(wQ)):
  731. for j in range(len(wQ)):
  732. term_QQ += wQ[i] * wQ[j] * k_gauss_of_gaussians(MUQ[i], SIGQ[i], MUQ[j], SIGQ[j], sigma_kernel)
  734. term_PQ = 0.0
  735. for i in range(len(wP)):
  736. for j in range(len(wQ)):
  737. term_PQ += wP[i] * wQ[j] * k_gauss_of_gaussians(MUP[i], SIGP[i], MUQ[j], SIGQ[j], sigma_kernel)
  739. dist2 = term_PP + term_QQ - 2.0 * term_PQ
  740. return max(float(dist2), 0.0)
  742. def compute_distance_matrix(models, sigma_kernel):
  743. n = len(models)
  744. D = np.zeros((n, n), dtype=float)
  746. print(f"Computing {n}x{n} RKHS distance matrix (σ_kernel={sigma_kernel:.6g})...")
  747. comps = [model_to_components(m) for m in models]
  749. def dist2_from_comps(cP, cQ):
  750. wP, MUP, SIGP = cP
  751. wQ, MUQ, SIGQ = cQ
  752. term_PP = 0.0
  753. for i in range(len(wP)):
  754. for j in range(len(wP)):
  755. term_PP += wP[i] * wP[j] * k_gauss_of_gaussians(MUP[i], SIGP[i], MUP[j], SIGP[j], sigma_kernel)
  756. term_QQ = 0.0
  757. for i in range(len(wQ)):
  758. for j in range(len(wQ)):
  759. term_QQ += wQ[i] * wQ[j] * k_gauss_of_gaussians(MUQ[i], SIGQ[i], MUQ[j], SIGQ[j], sigma_kernel)
  760. term_PQ = 0.0
  761. for i in range(len(wP)):
  762. for j in range(len(wQ)):
  763. term_PQ += wP[i] * wQ[j] * k_gauss_of_gaussians(MUP[i], SIGP[i], MUQ[j], SIGQ[j], sigma_kernel)
  764. return max(float(term_PP + term_QQ - 2.0 * term_PQ), 0.0)
  766. for i in range(n):
  767. D[i, i] = 0.0
  768. for j in range(i + 1, n):
  769. d2 = dist2_from_comps(comps[i], comps[j])
  770. D[i, j] = D[j, i] = d2
  771. if (i + 1) % 20 == 0:
  772. print(f" Processed {i + 1}/{n} rows")
  774. print(f"Done. Stats: min={D.min():.6f}, max={D.max():.6f}, mean={D.mean():.6f}")
  775. return D
  777. # Distance to kernel transformation
  778. def distance_to_kernel(D, gamma, epsilon=0.0):
  779. """
  780. Transform distance matrix to kernel (similarity) matrix using exponential transformation.
  782. Args:
  783. D: Distance matrix (n x n)
  784. gamma: Kernel bandwidth parameter
  785. epsilon: Optional diagonal regularization (default 0.0)
  787. Returns:
  788. K: Kernel matrix (n x n)
  789. """
  790. n = D.shape[0]
  791. K = np.exp(-D / (2 * gamma**2))
  793. if epsilon > 0:
  794. K += epsilon * np.eye(n)
  796. return K
  798. # Kernel matrix diagnostics
  799. def check_kernel_psd(K, verbose=True):
  800. """
  801. Check if kernel matrix is positive semi-definite.
  803. Args:
  804. K: Kernel matrix
  805. verbose: Print diagnostics
  807. Returns:
  808. is_psd: Boolean indicating if K is PSD
  809. min_eigval: Minimum eigenvalue
  810. """
  811. eigvals = np.linalg.eigvalsh(K)
  812. min_eigval = eigvals.min()
  813. max_eigval = eigvals.max()
  814. condition = max_eigval / max(abs(min_eigval), 1e-10)
  816. if verbose:
  817. print(f"Kernel diagnostics:")
  818. print(f" Min eigenvalue: {min_eigval:.6e}")
  819. print(f" Max eigenvalue: {max_eigval:.6e}")
  820. print(f" Condition number: {condition:.6e}")
  822. if min_eigval < -1e-10:
  823. print(f" WARNING: Kernel has negative eigenvalues!")
  824. elif min_eigval < 1e-12:
  825. print(f" WARNING: Very small positive eigenvalues, may be numerically unstable")
  826. else:
  827. print(f" Status: Kernel is well-conditioned")
  829. is_psd = min_eigval >= -1e-10
  830. return is_psd, min_eigval
  832. # ----------------------------
  833. # Cross-validation with SVM
  834. # ----------------------------
  835. def kfold_cross_validation_svm(models, labels, sigma_kernel, gamma, C,
  836. n_splits=5, random_state=42, epsilon=0.0):
  837. """
  838. K-fold cross-validation using SVM with precomputed kernel.
  840. Args:
  841. models: List of HMM-GMM models
  842. labels: Class labels
  843. sigma_kernel: Bandwidth for MMD² distance computation
  844. gamma: Bandwidth for kernel transformation
  845. C: SVM regularization parameter
  846. n_splits: Number of CV folds
  847. random_state: Random seed
  848. epsilon: Diagonal regularization (default 0.0)
  849. """
  850. skf = StratifiedKFold(n_splits=n_splits, shuffle=True, random_state=random_state)
  852. # Stage 1: Compute distance matrix
  853. full_D = compute_distance_matrix(models, sigma_kernel)
  855. # Stage 2: Transform to kernel matrix
  856. full_K = distance_to_kernel(full_D, gamma, epsilon)
  858. fold_results = {
  859. 'accuracy': [], 'balanced_accuracy': [], 'auc_roc': [], 'mcc': [],
  860. 'predictions': [], 'probabilities': [], 'true_labels': [],
  861. 'confusion_matrices': []
  862. }
  864. print(f"Performing {n_splits}-fold cross-validation with SVM...")
  865. print(f" σ_kernel={sigma_kernel:.6g}, γ={gamma:.6g}, C={C:.6g}")
  867. for fold, (train_idx, val_idx) in enumerate(skf.split(models, labels)):
  868. print(f" Fold {fold + 1}/{n_splits}")
  870. # Extract kernel submatrices for train and validation
  871. train_K = full_K[np.ix_(train_idx, train_idx)]
  872. val_train_K = full_K[np.ix_(val_idx, train_idx)]
  874. # Train SVM with precomputed kernel
  875. svm = SVC(kernel='precomputed', C=C, probability=True, random_state=random_state)
  876. svm.fit(train_K, labels[train_idx])
  878. # Predict on validation set
  879. val_pred = svm.predict(val_train_K)
  880. val_prob = svm.predict_proba(val_train_K)[:, 1]
  882. fold_results['predictions'].extend(val_pred.tolist())
  883. fold_results['probabilities'].extend(val_prob.tolist())
  884. fold_results['true_labels'].extend(labels[val_idx].tolist())
  886. # Compute metrics
  887. acc = accuracy_score(labels[val_idx], val_pred)
  888. bacc = balanced_accuracy_score(labels[val_idx], val_pred)
  889. auc_ = roc_auc_score(labels[val_idx], val_prob)
  890. mcc = matthews_corrcoef(labels[val_idx], val_pred)
  891. cm = confusion_matrix(labels[val_idx], val_pred)
  893. fold_results['accuracy'].append(acc)
  894. fold_results['balanced_accuracy'].append(bacc)
  895. fold_results['auc_roc'].append(auc_)
  896. fold_results['mcc'].append(mcc)
  897. fold_results['confusion_matrices'].append(cm)
  899. print(f" Acc={acc:.4f}, BalAcc={bacc:.4f}, AUC={auc_:.4f}, MCC={mcc:.4f}")
  901. # Convert to arrays
  902. for k in ['accuracy', 'balanced_accuracy', 'auc_roc', 'mcc']:
  903. fold_results[k] = np.array(fold_results[k])
  904. fold_results['predictions'] = np.array(fold_results['predictions'])
  905. fold_results['probabilities'] = np.array(fold_results['probabilities'])
  906. fold_results['true_labels'] = np.array(fold_results['true_labels'])
  908. return fold_results
  910. # -----------------------------
  911. # Optuna optimization for SVM
  912. # -----------------------------
  913. def hyperparameter_optimization_svm(models, labels,
  914. sigma_kernel_range=(1e-6, 1.0),
  915. gamma_range=(1e-6, 1.0),
  916. C_range=(1e-3, 1e3),
  917. n_trials=50, n_splits=5,
  918. timeout=1800, epsilon=0.0):
  919. """
  920. Hyperparameter optimization for SVM with RKHS kernel.
  922. Optimizes:
  923. - sigma_kernel: Bandwidth for MMD² distance computation
  924. - gamma: Bandwidth for kernel transformation
  925. - C: SVM regularization parameter
  926. """
  927. print("Starting hyperparameter optimization for SVM...")
  928. print(f"Search space:")
  929. print(f" sigma_kernel: log-uniform from {sigma_kernel_range[0]:.1e} to {sigma_kernel_range[1]:.1e}")
  930. print(f" gamma: log-uniform from {gamma_range[0]:.1e} to {gamma_range[1]:.1e}")
  931. print(f" C: log-uniform from {C_range[0]:.1e} to {C_range[1]:.1e}")
  932. print(f" n_trials: {n_trials}, CV folds: {n_splits}")
  934. def objective(trial):
  935. # Sample hyperparameters
  936. sigma_kernel = trial.suggest_float('sigma_kernel', sigma_kernel_range[0],
  937. sigma_kernel_range[1], log=True)
  938. gamma = trial.suggest_float('gamma', gamma_range[0], gamma_range[1], log=True)
  939. C = trial.suggest_float('C', C_range[0], C_range[1], log=True)
  941. try:
  942. # Stage 1: Compute distance matrix
  943. full_D = compute_distance_matrix(models, sigma_kernel)
  945. # Stage 2: Transform to kernel matrix
  946. full_K = distance_to_kernel(full_D, gamma, epsilon)
  948. # Cross-validation
  949. skf = StratifiedKFold(n_splits=n_splits, shuffle=True, random_state=42)
  950. scores = []
  952. for train_idx, val_idx in skf.split(models, labels):
  953. train_K = full_K[np.ix_(train_idx, train_idx)]
  954. val_train_K = full_K[np.ix_(val_idx, train_idx)]
  956. svm = SVC(kernel='precomputed', C=C, probability=True, random_state=42)
  957. svm.fit(train_K, labels[train_idx])
  959. pred = svm.predict(val_train_K)
  960. scores.append(balanced_accuracy_score(labels[val_idx], pred))
  962. return float(np.mean(scores))
  964. except Exception as e:
  965. print(f"Trial failed: {e}")
  966. return 0.0
  968. study = optuna.create_study(direction='maximize', sampler=TPESampler(seed=42),
  969. study_name="RKHS_SVM_Optimization")
  970. study.optimize(objective, n_trials=n_trials, timeout=timeout)
  972. print("\nOptimization completed!")
  973. print(f"Best balanced accuracy: {study.best_value:.4f}")
  974. print(f"Best parameters:")
  975. print(f" sigma_kernel: {study.best_params['sigma_kernel']:.6g}")
  976. print(f" gamma: {study.best_params['gamma']:.6g}")
  977. print(f" C: {study.best_params['C']:.6g}")
  979. return study
  981. # -----------------
  982. # Evaluation plots
  983. # -----------------
  984. def create_evaluation_plots(fold_results, class_names=['Control', 'ADHD']):
  985. fig, axes = plt.subplots(2, 3, figsize=(18, 12))
  987. overall_cm = confusion_matrix(fold_results['true_labels'], fold_results['predictions'])
  988. sns.heatmap(overall_cm, annot=True, fmt='d', cmap='Blues',
  989. xticklabels=class_names, yticklabels=class_names, ax=axes[0,0])
  990. axes[0,0].set_title('Overall Confusion Matrix')
  991. axes[0,0].set_xlabel('Predicted')
  992. axes[0,0].set_ylabel('Actual')
  994. fpr, tpr, _ = roc_curve(fold_results['true_labels'], fold_results['probabilities'])
  995. auc_score = auc(fpr, tpr)
  996. axes[0,1].plot(fpr, tpr, lw=2, label=f'ROC (AUC={auc_score:.3f})')
  997. axes[0,1].plot([0, 1], [0, 1], lw=2, linestyle='--', label='Random')
  998. axes[0,1].set_xlabel('FPR')
  999. axes[0,1].set_ylabel('TPR')
  1000. axes[0,1].set_title('ROC Curve')
  1001. axes[0,1].legend()
  1002. axes[0,1].grid(True, alpha=0.3)
  1004. precision, recall, _ = precision_recall_curve(fold_results['true_labels'],
  1005. fold_results['probabilities'])
  1006. pr_auc = auc(recall, precision)
  1007. axes[0,2].plot(recall, precision, lw=2, label=f'PR (AUC={pr_auc:.3f})')
  1008. axes[0,2].set_xlabel('Recall')
  1009. axes[0,2].set_ylabel('Precision')
  1010. axes[0,2].set_title('Precision-Recall')
  1011. axes[0,2].legend()
  1012. axes[0,2].grid(True, alpha=0.3)
  1014. metrics_data = [fold_results['accuracy'], fold_results['balanced_accuracy'],
  1015. fold_results['auc_roc'], fold_results['mcc']]
  1016. axes[1,0].boxplot(metrics_data, labels=['Acc', 'Bal.Acc', 'AUC', 'MCC'])
  1017. axes[1,0].set_title('CV Metric Distributions')
  1018. axes[1,0].set_ylabel('Score')
  1019. axes[1,0].grid(True, alpha=0.3)
  1021. folds = np.arange(1, len(fold_results['accuracy']) + 1)
  1022. axes[1,1].plot(folds, fold_results['balanced_accuracy'], 'o-', label='Balanced Acc')
  1023. axes[1,1].plot(folds, fold_results['auc_roc'], 'o-', label='AUC-ROC')
  1024. axes[1,1].set_xlabel('Fold')
  1025. axes[1,1].set_ylabel('Score')
  1026. axes[1,1].set_title('Performance by Fold')
  1027. axes[1,1].legend()
  1028. axes[1,1].grid(True, alpha=0.3)
  1030. class_counts = np.bincount(fold_results['true_labels'])
  1031. axes[1,2].pie(class_counts, labels=class_names, autopct='%1.1f%%', startangle=90)
  1032. axes[1,2].set_title('Class Distribution')
  1034. plt.tight_layout()
  1035. plt.show()
  1037. # -----------------------
  1038. # Full analysis pipeline with SVM
  1039. # -----------------------
  1040. def run_complete_analysis_svm(csv_path, test_size=0.2, random_state=42,
  1041. sigma_kernel_range=(1e-6, 1.0),
  1042. gamma_range=(1e-6, 1.0),
  1043. C_range=(1e-3, 1e3),
  1044. n_trials=50, n_splits=5, timeout=1800,
  1045. n_states=3, n_components=3, n_features=None,
  1046. epsilon=0.0, check_psd=True):
  1047. """
  1048. Complete analysis pipeline using SVM with RKHS kernel.
  1050. Args:
  1051. csv_path: Path to CSV with HMM-GMM parameters
  1052. test_size: Fraction for test set
  1053. random_state: Random seed
  1054. sigma_kernel_range: Range for distance computation bandwidth
  1055. gamma_range: Range for kernel transformation bandwidth
  1056. C_range: Range for SVM regularization
  1057. n_trials: Number of Optuna trials
  1058. n_splits: Number of CV folds
  1059. timeout: Optimization timeout in seconds
  1060. n_states: Number of HMM states
  1061. n_components: Number of GMM components per state
  1062. n_features: Number of features (auto-detected if None)
  1063. epsilon: Diagonal regularization for kernel matrix
  1064. check_psd: Whether to check kernel PSD property
  1065. """
  1066. print("=== RKHS-SVM ANALYSIS FOR HMM-GMM MODELS ===\n")
  1068. print("Loading CSV...")
  1069. df = pd.read_csv(csv_path)
  1070. df['class_label'] = (df['dataset_type'] == 'ADHD').astype(int)
  1072. print(f"Dataset summary:")
  1073. print(f" Total: {len(df)} | ADHD: {df['class_label'].sum()} | Control: {(1-df['class_label']).sum()}")
  1074. print(f" Class balance: {df['class_label'].mean():.3f}")
  1076. train_df, test_df = train_test_split(
  1077. df, test_size=test_size, random_state=random_state, stratify=df['class_label']
  1078. )
  1080. if n_features is None:
  1081. n_features = infer_n_features(train_df, n_states, n_components)
  1083. print("\nExtracting HMM-GMM parameters...")
  1084. train_models = extract_hmm_parameters_from_csv(train_df, n_states, n_components, n_features)
  1085. train_labels = train_df['class_label'].values
  1086. test_models = extract_hmm_parameters_from_csv(test_df, n_states, n_components, n_features)
  1087. test_labels = test_df['class_label'].values
  1089. # Hyperparameter optimization
  1090. study = hyperparameter_optimization_svm(
  1091. train_models, train_labels,
  1092. sigma_kernel_range=sigma_kernel_range,
  1093. gamma_range=gamma_range,
  1094. C_range=C_range,
  1095. n_trials=n_trials, n_splits=n_splits, timeout=timeout,
  1096. epsilon=epsilon
  1097. )
  1098. best_params = study.best_params
  1100. # Cross-validation with best parameters
  1101. print("\nCross-validation with best params...")
  1102. fold_results = kfold_cross_validation_svm(
  1103. train_models, train_labels,
  1104. best_params['sigma_kernel'],
  1105. best_params['gamma'],
  1106. best_params['C'],
  1107. n_splits=n_splits,
  1108. random_state=random_state,
  1109. epsilon=epsilon
  1110. )
  1112. print("\n" + "="*60)
  1113. print("CROSS-VALIDATION RESULTS")
  1114. print("="*60)
  1115. print(f" σ_kernel: {best_params['sigma_kernel']:.6g}")
  1116. print(f" γ (gamma): {best_params['gamma']:.6g}")
  1117. print(f" C: {best_params['C']:.6g}")
  1118. for name, key in [('Accuracy', 'accuracy'), ('Balanced Accuracy', 'balanced_accuracy'),
  1119. ('AUC-ROC','auc_roc'), ('MCC','mcc')]:
  1120. vals = fold_results[key]
  1121. print(f" {name}: {vals.mean():.4f} ± {vals.std():.4f}")
  1123. # Final test set evaluation
  1124. print("\n" + "="*60)
  1125. print("FINAL TEST SET EVALUATION")
  1126. print("="*60)
  1128. # Compute distance and kernel matrices
  1129. train_D = compute_distance_matrix(train_models, best_params['sigma_kernel'])
  1130. train_K = distance_to_kernel(train_D, best_params['gamma'], epsilon)
  1132. if check_psd:
  1133. print("\nChecking training kernel matrix...")
  1134. check_kernel_psd(train_K)
  1136. all_models = train_models + test_models
  1137. full_D = compute_distance_matrix(all_models, best_params['sigma_kernel'])
  1138. full_K = distance_to_kernel(full_D, best_params['gamma'], epsilon)
  1140. n_train = len(train_models)
  1141. test_train_K = full_K[n_train:, :n_train]
  1143. # Train final SVM on full training set
  1144. final_svm = SVC(kernel='precomputed', C=best_params['C'],
  1145. probability=True, random_state=random_state)
  1146. final_svm.fit(train_K, train_labels)
  1148. # Predict on test set
  1149. test_pred = final_svm.predict(test_train_K)
  1150. test_prob = final_svm.predict_proba(test_train_K)[:, 1]
  1152. # Compute test metrics
  1153. test_accuracy = accuracy_score(test_labels, test_pred)
  1154. test_bal_acc = balanced_accuracy_score(test_labels, test_pred)
  1155. test_auc = roc_auc_score(test_labels, test_prob)
  1156. test_mcc = matthews_corrcoef(test_labels, test_pred)
  1157. test_cm = confusion_matrix(test_labels, test_pred)
  1159. print(f" Accuracy: {test_accuracy:.4f}")
  1160. print(f" Balanced Accuracy: {test_bal_acc:.4f}")
  1161. print(f" AUC-ROC: {test_auc:.4f}")
  1162. print(f" Matthews Corr: {test_mcc:.4f}")
  1163. print("\nTest Confusion Matrix:")
  1164. print(f" Control ADHD")
  1165. print(f"Control {test_cm[0,0]:<7} {test_cm[0,1]:<7}")
  1166. print(f"ADHD {test_cm[1,0]:<7} {test_cm[1,1]:<7}")
  1168. # Create evaluation plots
  1169. create_evaluation_plots(fold_results)
  1171. return {
  1172. 'study': study,
  1173. 'best_params': best_params,
  1174. 'fold_results': fold_results,
  1175. 'test_results': {
  1176. 'accuracy': test_accuracy,
  1177. 'balanced_accuracy': test_bal_acc,
  1178. 'auc_roc': test_auc,
  1179. 'mcc': test_mcc,
  1180. 'confusion_matrix': test_cm,
  1181. 'predictions': test_pred,
  1182. 'probabilities': test_prob
  1183. },
  1184. 'final_model': final_svm
  1185. }
  1187. if __name__ == "__main__":
  1188. results = run_complete_analysis_svm(
  1189. # csv_path='.../synthetic_hmms_results.csv',
  1190. csv_path='.../NxGx/hmm_results.csv',
  1192. test_size=0.2,
  1193. random_state=42,
  1195. # Hyperparameter ranges
  1196. sigma_kernel_range=(1e-4, 1e-1), # For distance computation
  1197. gamma_range=(1e-5, 1e-2), # For kernel transformation
  1198. C_range=(1e-3, 1e+3), # SVM regularization
  1200. # Optimization settings
  1201. n_trials=250,
  1202. n_splits=7, # k-fold
  1203. timeout=(60 * 60 * 4), # hours
  1205. # Model configuration
  1206. n_states=3,# 4, 5
  1207. n_components=3, # 4, 5
  1208. n_features=None, # Auto-detect
  1210. # Kernel settings
  1211. epsilon=0.0, # No regularization by default
  1212. check_psd=True # Check kernel properties
  1213. )
  1215. print("\n" + "="*60)
  1216. print("ANALYSIS COMPLETE")
  1217. print("="*60)
  1218. print(f"Best σ_kernel: {results['best_params']['sigma_kernel']:.6g}")
  1219. print(f"Best γ (gamma): {results['best_params']['gamma']:.6g}")
  1220. print(f"Best C: {results['best_params']['C']:.6g}")
  1221. print(f"Final test balanced accuracy: {results['test_results']['balanced_accuracy']:.4f}")

rkhs-knn-svm-eeg-dataset.ipynb at commit d2eddf4, under MIT · at the source

Overview

Authors: Leonardo Lopez-Ortiz1, Cristhian K Valencia-Marin1, Julián Gil-González1, Paula M Herrera-Gómez2, David Cárdenas-Peña1
ORCID iDs: Cristhian K Valencia-Marin, Julián Gil-González, Paula M Herrera-Gómez, David Cárdenas-Peña
  1. Automatics Research Group, Universidad Tecnológica de Pereira (UTP), Pereira 660003, Colombia; (C.K.V.-M.); (J.G.-G.)
  2. Psiquiatría Neurociencias y Comunidad, Universidad Tecnológica de Pereira (UTP), Pereira 660003, Colombia
Institutions: Universidad Tecnológica de Pereira (Colombia)
Journal: Sensors (Basel, Switzerland), volume 26, issue 15, article 4773
Dates: received 19 May 2026; accepted 22 July 2026; published online 27 July 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.3390/s26154773 · PMID 42590549 · PMCID PMC13469683 · OpenAlex W7171436422
Open access: gold, a free copy (OpenAlex)
Status: code verified (how the code is found and verified)
Categories: EEG (modality), human (organism), ADHD (population)
Methods: Spectral & time-frequency, Preprocessing, Machine learning, Statistics, Physiology & signal measures, Smoothing, state filtering, decompositions
Keywords: EEG, ADHD, hidden Markov model, Gaussian mixture model, RKHS embedding, kernel classification
MeSH: Attention Deficit Disorder with Hyperactivity*, Electroencephalography*, Algorithms, Hidden Markov Models, Humans, Markov Chains (* major topic)
Topic: Attention Deficit Hyperactivity Disorder (Psychiatry and Mental health, Medicine), according to OpenAlex
Funding: Ministerio de Ciencia, Tecnología e Innovación (111091991908)
Citations: not cited yet (Europe PMC); 42 references in the paper

Abstract

Electroencephalography (EEG) provides a non-invasive means of supporting Attention-Deficit/Hyperactivity Disorder (ADHD) assessment. Nonetheless, existing pipelines often rely on handcrafted descriptors, segment-wise decisions, or deep architectures with limited subject-level generalization. This work introduces Hidden Markov Model-Induced Stationary RKHS Distance Learning (HIS), a probabilistic framework that represents each subject by a Hidden Markov Model with Gaussian-mixture emissions trained directly from frontal EEG recordings. Rather than vectorizing model parameters, each HMM is mapped to its induced stationary observation distribution and embedded into a Reproducing Kernel Hilbert Space (RKHS), where pairwise subject similarities are computed through a closed-form Hilbert embedding distance. These similarities are subsequently exploited by precomputed-kernel classifiers for subject-level prediction. The proposed method was evaluated against the Probability Product Kernel baseline using both a controlled synthetic EEG benchmark and a public pediatric ADHD dataset under progressively more rigorous validation protocols, culminating in repeated nested cross-validation with bootstrap confidence intervals and permutation testing. On the synthetic benchmark, HIS achieved 95.0% held-out accuracy and consistently outperformed the baseline across classifiers. On a real EEG dataset with 121 subjects, the primary evaluation protocol yielded a balanced accuracy of 73.5% (95% CI: 69.8–77.0%), an AUC of 79.6%, and an MCC of 0.483 (permutation p < 0.001) using an SVM with compact subject-specific HMMs. Complementary hyperparameter analyses and t-SNE visualizations demonstrated that HIS induces more stable and discriminative subject representations than the baseline. These results establish stationary RKHS embeddings of subject-specific HMMs as a leakage-aware framework for EEG-based ADHD decision support and underscore the critical influence of statistically rigorous evaluation protocols on reported classification performance.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 12 matches between paragraphs and lines of code.

leonlpz/HMM_RKHS_Classification_Project

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: d2eddf4778b0944543c6ea141c817cdd2323cc00, 11 July 2026
Languages: Jupyter (6), Python (1)
Size: 1,158 files, 7 scripts
Software Heritage: not archived
Found in: “Data Availability Statement”
Holds: README, license file, 6 notebooks
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (7 files), pandas (7 files), scikit-learn (7 files), Matplotlib (5 files), SciPy (5 files), seaborn (4 files), Numba (3 files), TensorFlow (3 files)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
9 files
At the source: github.com/leonlpz/HMM_RKHS_Classification_Project

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 7 scripts, each with its path and the digest of its content;
  • 12 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data Availability Statement

The dataset EEG data for ADHD / Control children, used in the validation of the current study, is publicly available at https://ieee-dataport.org/open-access/eeg-data-adhd-control-children. The code, baseline reimplementation, full per-configuration results, and trained model artifacts supporting the findings of this study are openly available in the project’s repository (https://github.com/leonlpz/HMM_RKHS_Classification_Project (accessed on 19 May 2026)).

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 5 authors, 6 keywords, 6 MeSH terms, 1 funder, 39 references.

Cite

This paper

Lopez-Ortiz, L., Valencia-Marin, C. K., Gil-González, J., Herrera-Gómez, P. M., & Cárdenas-Peña, D. (2026). EEG-Based Supported Diagnosis of ADHD Using Subject-Specific HMMs and Stationary RKHS Embeddings. Sensors (Basel, Switzerland), 26(15), 4773. https://doi.org/10.3390/s26154773

BibTeX

@article{lopezortiz2026eeg,
author = {Lopez-Ortiz, Leonardo and Valencia-Marin, Cristhian K and Gil-González, Julián and Herrera-Gómez, Paula M and Cárdenas-Peña, David},
title = {{EEG-Based Supported Diagnosis of ADHD Using Subject-Specific HMMs and Stationary RKHS Embeddings}},
journal = {Sensors (Basel, Switzerland)},
year = {2026},
month = jul,
volume = {26},
number = {15},
pages = {4773},
publisher = {Multidisciplinary Digital Publishing Institute (MDPI)},
issn = {1424-8220},
doi = {10.3390/s26154773},
url = {https://doi.org/10.3390/s26154773},
pmid = {42590549},
pmcid = {PMC13469683}
}

RIS

TY - JOUR
AU - Lopez-Ortiz, Leonardo
AU - Valencia-Marin, Cristhian K
AU - Gil-González, Julián
AU - Herrera-Gómez, Paula M
AU - Cárdenas-Peña, David
TI - EEG-Based Supported Diagnosis of ADHD Using Subject-Specific HMMs and Stationary RKHS Embeddings
T2 - Sensors (Basel, Switzerland)
J2 - Sensors (Basel)
PY - 2026
DA - 2026/07/27
VL - 26
IS - 15
SP - 4773
SN - 1424-8220
PB - Multidisciplinary Digital Publishing Institute (MDPI)
DO - 10.3390/s26154773
UR - https://doi.org/10.3390/s26154773
LA - en
ER -

CSL-JSON

{
"id": "10.3390/s26154773",
"type": "article-journal",
"title": "EEG-Based Supported Diagnosis of ADHD Using Subject-Specific HMMs and Stationary RKHS Embeddings",
"container-title": "Sensors (Basel, Switzerland)",
"author": [
{
"family": "Lopez-Ortiz",
"given": "Leonardo"
},
{
"family": "Valencia-Marin",
"given": "Cristhian K"
},
{
"family": "Gil-González",
"given": "Julián"
},
{
"family": "Herrera-Gómez",
"given": "Paula M"
},
{
"family": "Cárdenas-Peña",
"given": "David"
}
],
"container-title-short": "Sensors (Basel)",
"volume": "26",
"issue": "15",
"page": "4773",
"DOI": "10.3390/s26154773",
"PMID": "42590549",
"PMCID": "PMC13469683",
"ISSN": "1424-8220",
"publisher": "Multidisciplinary Digital Publishing Institute (MDPI)",
"URL": "https://doi.org/10.3390/s26154773",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
27
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.3390/s23073763 [code]
A Novel OpenBCI Framework for EEG-Based Neurophysiological Experiments.
Journal: Sensors (Basel, Switzerland)
In common: TensorFlow, seaborn, scikit-learn, 3 other tools, EEG, 1 reference, author David Cárdenas-Peña
[2] doi:10.1016/j.heliyon.2024.e26028 [code]
A novel approach to identify the brain regions that best classify ADHD by means of EEG and deep learning.
Journal: Heliyon
In common: TensorFlow, scikit-learn, pandas, 3 other tools, ADHD, EEG, 3 references
[3] doi:10.1371/journal.pone.0357213
Integrating socioeconomic context with multimodal EEG data for improved ADHD risk screening.
Journal: PloS one
In common: ADHD, EEG, 5 references
[4] doi:10.3390/s25227067
Takens-Based Kernel Transfer Entropy Connectivity Network for Motor Imagery Classification.
Journal: Sensors (Basel, Switzerland)
In common: EEG, 2 references, author David Cárdenas-Peña
[5] doi:10.21203/rs.3.rs-9010375/v1 [code]
Interpretable metrics for evaluating fidelity, diversity, and privacy in synthetic EEG data
Journal: Research Square (preprint)
In common: Numba, seaborn, scikit-learn, 4 other tools, EEG, 1 reference
[6] doi:10.1371/journal.pone.0347671 [code]
RMETNet: A cross-subject motor imagery EEG signal classification model based on TSLANet and riemannian geometry features.
Journal: PloS one
In common: TensorFlow, seaborn, scikit-learn, 4 other tools, EEG, 1 reference
[7] doi:10.1038/s41598-024-62230-9 [code]
Feasibility of decoding covert speech in ECoG with a Transformer trained on overt speech.
Journal: Scientific reports
In common: TensorFlow, seaborn, scikit-learn, 4 other tools, 1 reference
[8] doi:10.1007/s10484-024-09675-w [code]
Theta/Beta Ratio Neurofeedback Effects on Resting and Task-Related Theta Activity in Children with ADHD.
Journal: Applied psychophysiology and biofeedback
In common: seaborn, pandas, SciPy, 2 other tools, ADHD, EEG, 1 reference
[9] doi:10.1038/s41598-026-42310-8 [code]
Sleep awake detection from leg-worn wearables using deep sensor fusion.
Journal: Scientific reports
In common: TensorFlow, seaborn, scikit-learn, 4 other tools, ADHD
[10] doi:10.1038/s41598-026-50791-w [code]
Hybrid Vi+ECNN framework for advanced ADHD diagnostic accuracy in medical imaging.
Journal: Scientific reports
In common: seaborn, scikit-learn, pandas, 2 other tools, ADHD, 1 reference

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Add the badge to the code's README

OSCR: the badge of this paper

If you are an author or a maintainer of this paper's code, you can add this badge. It says what OSCR knows of the paper (its code, its validated tracing map) and links the code to this page. Paste one of these snippets under the README's title. Keep the two marks around it, so that the badge can be replaced in place and is never added twice. OSCR changes nothing in your repository.

Markdown
<!-- oscr-badge:start -->
[![OSCR](https://openscicode.org/badge/doi/10.3390/s26154773.svg)](https://openscicode.org/badge/doi/10.3390/s26154773)
<!-- oscr-badge:end -->
reStructuredText
.. oscr-badge:start

.. image:: https://openscicode.org/badge/doi/10.3390/s26154773.svg
   :target: https://openscicode.org/badge/doi/10.3390/s26154773
   :alt: OSCR

.. oscr-badge:end
HTML
<!-- oscr-badge:start -->
<a href="https://openscicode.org/badge/doi/10.3390/s26154773"><img src="https://openscicode.org/badge/doi/10.3390/s26154773.svg" alt="OSCR"></a>
<!-- oscr-badge:end -->

On GitHub, open the README, choose Edit (the pencil), paste the snippet under the title, then "Commit changes".

If you choose to, you can have the badge kept up to date for you, with oscr badge --write from the command line or with the GitHub Action in your repository's workflows. Every other form and option (a card, AsciiDoc, LaTeX, shields.io) is on the badge page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: code errors and mismatches between the code and the paper about this paper and its code, from readers and its authors, beside the paper and its tracing map.

Report a mismatch between the code and the paper · Report a code error (signed in)

Reproductions: reports from readers who ran the authors' code and did not get the paper's result, what they ran, with which environment, commit and data.

Report a reproduction failure (signed in)

Versions of its code: the releases of its code tied to a version of this paper (the preprint, the accepted manuscript, the version of record…), each with the tracing map versioned with it.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.