Refactored pyvbmc code

This commit is contained in:
Louis Lacoste 2025-12-07 17:38:03 +01:00
parent ba4e88d54a
commit a3e9eb68eb

View file

@ -1,3 +1,4 @@
# %%
# ============================================================================
# OPTIMISATION PYVBMC POUR MODÈLES Q-LEARNING AVEC ÉVÉNEMENTS RARES
# ============================================================================
@ -14,6 +15,8 @@ from pathlib import Path
# Tentative d'import PyVBMC
try:
from pyvbmc import VBMC
from pyvbmc.priors import UniformBox
PYVBMC_AVAILABLE = True
except ImportError:
PYVBMC_AVAILABLE = False
@ -26,6 +29,8 @@ from load_data import all_participant_data, unique_participants
# CONFIGURATIONS DES MODÈLES EMBOÎTÉS
# ============================================================================
# %%
def get_model_configs() -> Dict:
"""Retourne les configurations des différents modèles."""
return {
@ -38,7 +43,7 @@ def get_model_configs() -> Dict:
"n_params": 3,
"param_names": ["alpha", "forget", "lambda"],
"lower": np.array([-5, -5, -3]),
"upper": np.array([5, 5, 3])
"upper": np.array([5, 5, 3]),
},
"GAIN_LOSS": {
"name": "GAIN_LOSS",
@ -49,7 +54,7 @@ def get_model_configs() -> Dict:
"n_params": 4,
"param_names": ["alpha_loss", "alpha_gain", "forget", "lambda"],
"lower": np.array([-5, -5, -5, -3]),
"upper": np.array([5, 5, 5, 3])
"upper": np.array([5, 5, 5, 3]),
},
"BIASED": {
"name": "BIASED",
@ -59,12 +64,19 @@ def get_model_configs() -> Dict:
"has_rho": False,
"n_params": 10,
"param_names": [
"alpha_loss", "alpha_gain",
"forget_1", "forget_2", "forget_3", "forget_4",
"lambda_1", "lambda_2", "lambda_3", "lambda_4"
"alpha_loss",
"alpha_gain",
"forget_1",
"forget_2",
"forget_3",
"forget_4",
"lambda_1",
"lambda_2",
"lambda_3",
"lambda_4",
],
"lower": np.concatenate([[-5, -5], np.full(4, -5), np.full(4, -3)]),
"upper": np.concatenate([[5, 5], np.full(4, 5), np.full(4, 3)])
"upper": np.concatenate([[5, 5], np.full(4, 5), np.full(4, 3)]),
},
"REE_BIASED_SIMPLE": {
"name": "REE_BIASED_SIMPLE",
@ -74,11 +86,15 @@ def get_model_configs() -> Dict:
"has_rho": True,
"n_params": 6,
"param_names": [
"alpha_loss", "alpha_gain", "forget", "lambda",
"rho_BS", "rho_JP"
"alpha_loss",
"alpha_gain",
"forget",
"lambda",
"rho_BS",
"rho_JP",
],
"lower": np.array([-5, -5, -5, -3, -10, -10]),
"upper": np.array([5, 5, 5, 3, 10, 10])
"upper": np.array([5, 5, 5, 3, 10, 10]),
},
"REE_BIASED_COMPLEX": {
"name": "REE_BIASED_COMPLEX",
@ -88,13 +104,23 @@ def get_model_configs() -> Dict:
"has_rho": True,
"n_params": 12,
"param_names": [
"alpha_loss", "alpha_gain",
"forget_1", "forget_2", "forget_3", "forget_4",
"lambda_1", "lambda_2", "lambda_3", "lambda_4",
"rho_BS", "rho_JP"
"alpha_loss",
"alpha_gain",
"forget_1",
"forget_2",
"forget_3",
"forget_4",
"lambda_1",
"lambda_2",
"lambda_3",
"lambda_4",
"rho_BS",
"rho_JP",
],
"lower": np.concatenate([[-5, -5], np.full(4, -5), np.full(4, -3), [-10, -10]]),
"upper": np.concatenate([[5, 5], np.full(4, 5), np.full(4, 3), [10, 10]])
"lower": np.concatenate(
[[-5, -5], np.full(4, -5), np.full(4, -3), [-10, -10]]
),
"upper": np.concatenate([[5, 5], np.full(4, 5), np.full(4, 3), [10, 10]]),
},
"REE_LEARNING_SIMPLE": {
"name": "REE_LEARNING_SIMPLE",
@ -104,11 +130,15 @@ def get_model_configs() -> Dict:
"has_rho": False,
"n_params": 6,
"param_names": [
"alpha_loss", "alpha_gain", "alpha_BS", "alpha_JP",
"forget", "lambda"
"alpha_loss",
"alpha_gain",
"alpha_BS",
"alpha_JP",
"forget",
"lambda",
],
"lower": np.array([-5, -5, -5, -5, -5, -3]),
"upper": np.array([5, 5, 5, 5, 5, 3])
"upper": np.array([5, 5, 5, 5, 5, 3]),
},
"REE_LEARNING_COMPLEX": {
"name": "REE_LEARNING_COMPLEX",
@ -118,12 +148,21 @@ def get_model_configs() -> Dict:
"has_rho": False,
"n_params": 12,
"param_names": [
"alpha_loss", "alpha_gain", "alpha_BS", "alpha_JP",
"forget_1", "forget_2", "forget_3", "forget_4",
"lambda_1", "lambda_2", "lambda_3", "lambda_4"
"alpha_loss",
"alpha_gain",
"alpha_BS",
"alpha_JP",
"forget_1",
"forget_2",
"forget_3",
"forget_4",
"lambda_1",
"lambda_2",
"lambda_3",
"lambda_4",
],
"lower": np.concatenate([[-5, -5, -5, -5], np.full(4, -5), np.full(4, -3)]),
"upper": np.concatenate([[5, 5, 5, 5], np.full(4, 5), np.full(4, 3)])
"upper": np.concatenate([[5, 5, 5, 5], np.full(4, 5), np.full(4, 3)]),
},
"REE_LEARNING_BIASED_SIMPLE": {
"name": "REE_LEARNING_BIASED_SIMPLE",
@ -133,11 +172,17 @@ def get_model_configs() -> Dict:
"has_rho": True,
"n_params": 8,
"param_names": [
"alpha_loss", "alpha_gain", "alpha_BS", "alpha_JP",
"forget", "lambda", "rho_BS", "rho_JP"
"alpha_loss",
"alpha_gain",
"alpha_BS",
"alpha_JP",
"forget",
"lambda",
"rho_BS",
"rho_JP",
],
"lower": np.array([-5, -5, -5, -5, -5, -3, -10, -10]),
"upper": np.array([5, 5, 5, 5, 5, 3, 10, 10])
"upper": np.array([5, 5, 5, 5, 5, 3, 10, 10]),
},
"REE_LEARNING_BIASED_COMPLEX": {
"name": "REE_LEARNING_BIASED_COMPLEX",
@ -147,14 +192,28 @@ def get_model_configs() -> Dict:
"has_rho": True,
"n_params": 14,
"param_names": [
"alpha_loss", "alpha_gain", "alpha_BS", "alpha_JP",
"forget_1", "forget_2", "forget_3", "forget_4",
"lambda_1", "lambda_2", "lambda_3", "lambda_4",
"rho_BS", "rho_JP"
"alpha_loss",
"alpha_gain",
"alpha_BS",
"alpha_JP",
"forget_1",
"forget_2",
"forget_3",
"forget_4",
"lambda_1",
"lambda_2",
"lambda_3",
"lambda_4",
"rho_BS",
"rho_JP",
],
"lower": np.concatenate([[-5, -5, -5, -5], np.full(4, -5), np.full(4, -3), [-10, -10]]),
"upper": np.concatenate([[5, 5, 5, 5], np.full(4, 5), np.full(4, 3), [10, 10]])
}
"lower": np.concatenate(
[[-5, -5, -5, -5], np.full(4, -5), np.full(4, -3), [-10, -10]]
),
"upper": np.concatenate(
[[5, 5, 5, 5], np.full(4, 5), np.full(4, 3), [10, 10]]
),
},
}
@ -162,8 +221,13 @@ def get_model_configs() -> Dict:
# MODÈLE Q-LEARNING GÉNÉRIQUE
# ============================================================================
def qlearning_generic(params: np.ndarray, data: pd.DataFrame, model_config: Dict,
return_negLL: bool = True) -> float:
def qlearning_generic(
params: np.ndarray,
data: pd.DataFrame,
model_config: Dict,
return_negLL: bool = True,
) -> float:
"""
Modèle Q-learning générique avec support pour différentes architectures de paramètres.
@ -204,7 +268,7 @@ def qlearning_generic(params: np.ndarray, data: pd.DataFrame, model_config: Dict
forget = np.full(n_arms, expit(params[param_idx]))
param_idx += 1
elif model_config["n_forget"] == 4:
forget = expit(params[param_idx:(param_idx + 4)])
forget = expit(params[param_idx : (param_idx + 4)])
param_idx += 4
# LAMBDA(S)
@ -212,7 +276,7 @@ def qlearning_generic(params: np.ndarray, data: pd.DataFrame, model_config: Dict
lambda_vals = np.full(n_arms, np.exp(params[param_idx]))
param_idx += 1
elif model_config["n_lambda"] == 4:
lambda_vals = np.exp(params[param_idx:(param_idx + 4)])
lambda_vals = np.exp(params[param_idx : (param_idx + 4)])
param_idx += 4
# RHO(S) - Biais pour événements rares
@ -277,12 +341,15 @@ def qlearning_generic(params: np.ndarray, data: pd.DataFrame, model_config: Dict
return log_lik
# %%
# ============================================================================
# OPTIMISATION AVEC PYVBMC
# ============================================================================
def fit_participant_pyvbmc(participant_data: pd.DataFrame, model_config: Dict,
verbose: bool = True) -> Dict:
def fit_participant_pyvbmc(
participant_data: pd.DataFrame, model_config: Dict, verbose: bool = True
) -> Dict:
"""
Optimise les paramètres du modèle pour un participant utilisant PyVBMC.
@ -295,13 +362,17 @@ def fit_participant_pyvbmc(participant_data: pd.DataFrame, model_config: Dict,
Dictionnaire avec les résultats d'optimisation
"""
if not PYVBMC_AVAILABLE:
raise RuntimeError("PyVBMC n'est pas installé. Installez avec: pip install pyvbmc")
raise RuntimeError(
"PyVBMC n'est pas installé. Installez avec: pip install pyvbmc"
)
# Définition de la fonction de log-densité pour PyVBMC
def log_posterior(params_array):
def log_likelihood(params_array):
"""PyVBMC maximise, donc on retourne -negLL."""
params = np.asarray(params_array).flatten()
negLL = qlearning_generic(params, participant_data, model_config, return_negLL=True)
negLL = qlearning_generic(
params, participant_data, model_config, return_negLL=True
)
return -negLL
# Point de départ (milieu des bornes)
@ -319,22 +390,25 @@ def fit_participant_pyvbmc(participant_data: pd.DataFrame, model_config: Dict,
# Initialisation et optimisation de VBMC
vbmc = VBMC(
log_posterior,
log_likelihood,
x0,
model_config["lower"],
model_config["upper"],
plb,
pub,
options={
"verbose": 0 if not verbose else 1,
# "verbose": 0 if not verbose else 1,
"display": "off",
}
},
prior=UniformBox(
a=model_config["lower"], b=model_config["upper"], D=model_config["n_params"]
),
)
vp, results = vbmc.optimize()
# Extraction des statistiques
posterior_mean, posterior_cov = vp.moments()
posterior_mean, posterior_cov = vp.moments(orig_flag=True, cov_flag=True)
posterior_mean = np.asarray(posterior_mean).flatten()
posterior_sd = np.sqrt(np.diag(posterior_cov))
@ -344,7 +418,9 @@ def fit_participant_pyvbmc(participant_data: pd.DataFrame, model_config: Dict,
n_iterations = results.get("iterations", np.nan)
# Calcul du negLL avec la posterior mean
negLL = qlearning_generic(posterior_mean, participant_data, model_config, return_negLL=True)
negLL = qlearning_generic(
posterior_mean, participant_data, model_config, return_negLL=True
)
n_obs = len(participant_data)
# Calcul des critères d'information
@ -365,7 +441,7 @@ def fit_participant_pyvbmc(participant_data: pd.DataFrame, model_config: Dict,
"posterior_mean": posterior_mean,
"posterior_sd": posterior_sd,
"vp": vp,
"results": results
"results": results,
}
# Ajout des paramètres estimés
@ -376,8 +452,13 @@ def fit_participant_pyvbmc(participant_data: pd.DataFrame, model_config: Dict,
return result
def fit_participant_deoptim(participant_data: pd.DataFrame, model_config: Dict,
n_runs: int = 5, verbose: bool = True) -> Dict:
def fit_participant_deoptim(
participant_data: pd.DataFrame,
model_config: Dict,
n_runs: int = 5,
verbose: bool = True,
n_workers: int = 1,
) -> Dict:
"""
Optimise les paramètres du modèle pour un participant utilisant minimisation scipy.
@ -403,16 +484,18 @@ def fit_participant_deoptim(participant_data: pd.DataFrame, model_config: Dict,
np.random.seed(1000 * hash(model_config["name"]) % (2**31) + run)
def objective(params):
return qlearning_generic(params, participant_data, model_config, return_negLL=True)
return qlearning_generic(
params, participant_data, model_config, return_negLL=True
)
result = differential_evolution(
objective,
bounds=list(zip(model_config["lower"], model_config["upper"])),
maxiter=200,
popsize=max(50, model_config["n_params"] * 10),
seed=1000 * hash(model_config["name"]) % (2**31) + run,
workers=1,
updating="deferred"
rng=1000 * hash(model_config["name"]) % (2**31) + run,
workers=n_workers,
updating="deferred",
)
all_negLLs.append(result.fun)
@ -458,9 +541,14 @@ def fit_participant_deoptim(participant_data: pd.DataFrame, model_config: Dict,
# OPTIMISATION POUR TOUS LES PARTICIPANTS ET MODÈLES
# ============================================================================
def fit_all_participants(data: pd.DataFrame, models_to_fit: Optional[List[str]] = None,
method: str = "VBMC", n_participants: Optional[int] = None,
verbose: bool = True) -> Dict[str, List[Dict]]:
def fit_all_participants(
data: pd.DataFrame,
models_to_fit: Optional[List[str]] = None,
method: str = "VBMC",
n_participants: Optional[int] = None,
verbose: bool = True,
) -> Dict[str, List[Dict]]:
"""
Ajuste tous les modèles pour tous les participants.
@ -499,10 +587,13 @@ def fit_all_participants(data: pd.DataFrame, models_to_fit: Optional[List[str]]
try:
if method == "VBMC":
result = fit_participant_pyvbmc(participant_data, model_config, verbose=False)
result = fit_participant_pyvbmc(
participant_data, model_config, verbose=False
)
else:
result = fit_participant_deoptim(participant_data, model_config,
n_runs=5, verbose=False)
result = fit_participant_deoptim(
participant_data, model_config, n_runs=5, verbose=False
)
result["participant"] = participant_id
model_results.append(result)
@ -523,6 +614,7 @@ def fit_all_participants(data: pd.DataFrame, models_to_fit: Optional[List[str]]
# COMPARAISON DES MODÈLES
# ============================================================================
def compare_models(all_results: Dict[str, List[Dict]]) -> Dict:
"""
Compare les modèles et sélectionne les meilleurs par participant.
@ -563,16 +655,20 @@ def compare_models(all_results: Dict[str, List[Dict]]) -> Dict:
all_results_list = []
for model_name, results in all_results.items():
for result in results:
all_results_list.append({
"participant": result["participant"],
"model": model_name,
"BIC": result["BIC"],
"AIC": result["AIC"],
"negLL": result["negLL"]
})
all_results_list.append(
{
"participant": result["participant"],
"model": model_name,
"BIC": result["BIC"],
"AIC": result["AIC"],
"negLL": result["negLL"],
}
)
all_results_df = pd.DataFrame(all_results_list)
best_per_participant = all_results_df.loc[all_results_df.groupby("participant")["BIC"].idxmin()]
best_per_participant = all_results_df.loc[
all_results_df.groupby("participant")["BIC"].idxmin()
]
print("\n=== BEST MODELS PER PARTICIPANT ===")
print(best_per_participant["model"].value_counts())
@ -580,7 +676,7 @@ def compare_models(all_results: Dict[str, List[Dict]]) -> Dict:
return {
"global_comparison": global_comparison_df,
"best_per_participant": best_per_participant,
"all_results": all_results
"all_results": all_results,
}
@ -588,7 +684,10 @@ def compare_models(all_results: Dict[str, List[Dict]]) -> Dict:
# SAUVEGARDE DES RÉSULTATS
# ============================================================================
def save_results(all_results: Dict[str, List[Dict]], output_dir: str = "results") -> None:
def save_results(
all_results: Dict[str, List[Dict]], output_dir: str = "results"
) -> None:
"""
Sauvegarde les résultats d'optimisation en CSV.
@ -603,21 +702,136 @@ def save_results(all_results: Dict[str, List[Dict]], output_dir: str = "results"
results_df = pd.DataFrame(results)
# Garder seulement les colonnes numériques pour le CSV
cols_to_keep = [col for col in results_df.columns
if col not in ["vp", "results", "posterior_mean", "posterior_sd"]]
cols_to_keep = [
col
for col in results_df.columns
if col not in ["vp", "results", "posterior_mean", "posterior_sd"]
]
results_df[cols_to_keep].to_csv(
output_path / f"results_{model_name}.csv",
index=False
output_path / f"results_{model_name}.csv", index=False
)
print(f"Saved: results_{model_name}.csv")
def fit_vbmc_and_diffEvol(
participant_data: pd.DataFrame,
model_config: Dict,
n_deoptim_runs: int = 5,
n_workers: int = 1,
verbose: bool = True,
) -> Tuple[Dict, Dict]:
"""
Ajuste un modèle à l'aide de PyVBMC et Differential Evolution pour comparaison.
Args:
participant_data: Données du participant
model_config: Configuration du modèle
n_deoptim_runs: Nombre de runs pour Differential Evolution
verbose: Affiche les progressions
Returns:
Tuple avec les résultats VBMC et Differential Evolution
"""
if verbose:
print(f" Fitting with VBMC")
vbmc_result = fit_participant_pyvbmc(
participant_data, model_config, verbose=verbose
)
if verbose:
print(f" Fitting with Differential Evolution")
deoptim_result = fit_participant_deoptim(
participant_data,
model_config,
n_runs=n_deoptim_runs,
n_workers=n_workers,
verbose=verbose,
)
return vbmc_result, deoptim_result
def fit_all_participants_both_methods(
data: pd.DataFrame,
models_to_fit: Optional[List[str]] = None,
n_participants: Optional[int] = None,
n_deoptim_runs: int = 5,
n_workers: int = 1,
verbose: bool = True,
) -> Dict[str, List[Dict]]:
"""
Ajuste tous les modèles pour tous les participants avec les deux méthodes.
Args:
data: DataFrame avec les données de tous les participants
models_to_fit: Liste des noms de modèles à ajuster (None = tous)
n_participants: Nombre de participants à traiter (None = tous)
n_deoptim_runs: Nombre de runs pour Differential Evolution
verbose: Affiche les progressions
Returns:
Dictionnaire avec les résultats par modèle et méthode
"""
model_configs = get_model_configs()
if models_to_fit is not None:
model_configs = {k: v for k, v in model_configs.items() if k in models_to_fit}
participants = data["participant"].unique()
if n_participants is not None:
participants = participants[:n_participants]
all_results = {}
for model_name, model_config in model_configs.items():
if verbose:
print(f"\n=== Fitting model: {model_name} ===")
model_results = []
for participant_id in participants:
if verbose:
print(f" Participant: {participant_id}")
participant_data = data[data["participant"] == participant_id].copy()
try:
vbmc_result, deoptim_result = fit_vbmc_and_diffEvol(
participant_data,
model_config,
n_deoptim_runs=n_deoptim_runs,
n_workers=n_workers,
verbose=False,
)
vbmc_result["participant"] = participant_id
deoptim_result["participant"] = participant_id
model_results.append(
{"VBMC": vbmc_result, "Differential_Evolution": deoptim_result}
)
if verbose:
print(
f" VBMC negLL: {vbmc_result['negLL']:.2f}, BIC: {vbmc_result['BIC']:.2f}"
)
print(
f" DE negLL: {deoptim_result['negLL']:.2f}, BIC: {deoptim_result['BIC']:.2f}"
)
except Exception as e:
print(f" ERROR: {str(e)}")
continue
all_results[model_name] = model_results
return all_results
# %%
# ============================================================================
# EXEMPLE D'UTILISATION
# ============================================================================
if __name__ == "__main__":
print("=== PyVBMC Optimization for Q-Learning Models ===\n")
print("=== Optimization for Q-Learning Models ===\n")
# Préparation des données
print("Loading data...")
@ -634,14 +848,32 @@ if __name__ == "__main__":
print(f" PyVBMC not available - using {method}")
# Ajustement de quelques modèles pour test
models_to_fit = ["HOMOGENEOUS", "GAIN_LOSS", "REE_BIASED_SIMPLE"]
models_to_fit = [
"HOMOGENEOUS",
"GAIN_LOSS",
"REE_BIASED_SIMPLE",
"REE_BIASED_COMPLEX",
"REE_LEARNING_SIMPLE",
"REE_LEARNING_COMPLEX",
"REE_LEARNING_BIASED_SIMPLE",
"REE_LEARNING_BIASED_COMPLEX",
]
# all_results = fit_all_participants_both_methods(
# data_for_fitting,
# models_to_fit=models_to_fit,
# # method=method,
# n_participants=2, # Set to a number to limit for testing
# n_workers=1,
# verbose=True,
# )
all_results = fit_all_participants(
data_for_fitting,
models_to_fit=models_to_fit,
method=method,
n_participants=2, # Set to a number to limit for testing
verbose=True
n_participants=1, # Set to a number to limit for testing
verbose=True,
)
# Comparaison des modèles
@ -654,3 +886,5 @@ if __name__ == "__main__":
comparison["best_per_participant"].to_csv("results/best_models.csv", index=False)
print("\nDone!")
# %%