Exploring hierarchical estimation of hyperparameters (cpm.hierarchical)¶

import warnings
import ipyparallel as ipp  # for parallel computing with ipython (specific for Jupyter Notebook)

import cpm
import cpm.datasets as datasets
from cpm.generators import Parameters, Value

import numpy
import pandas as pd
import functions as f ## I prepped some models here, so you can use them directly

## plotting libraries
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import seaborn as sns

%matplotlib inline
## Hierarchical fitting can take a lot of time depending on your setup, so I ran the simulation for you already, so we can focus on exploring the results. 
hyperparameters = pd.read_csv('04-fitting-hierarchical-hyperparameters.csv')
hyperparameters.head()

plt.figure(figsize=(10, 6))
sns.lineplot(
    data=hyperparameters.drop_duplicates(subset=['chain', 'iteration']),
    x='iteration',
    y='lme',
    hue='chain',
    marker='o',
    markersize=10,
    palette='tab10',
)
plt.title('Log Model Evidence (lme) over Iterations for Each Chain')
plt.xlabel('Iteration')
plt.ylabel('Log Model Evidence (lme)')
plt.legend(title='Chain')
plt.tight_layout()

fig, axes = plt.subplots(2, 1, figsize=(7, 10), sharey=False)

sns.lineplot(
    data=hyperparameters[hyperparameters['parameter'] == 'alpha'],
    x='iteration',
    y='mean',
    hue='chain',
    marker='o',
    palette='tab10',
    ax=axes[0]
)
axes[0].set_title('Alpha Hyperparameter Mean by Iteration')
axes[0].set_xlabel('Iteration')
axes[0].set_ylabel('Mean (alpha)')
axes[0].legend(title='Chain')
axes[0].set_ylim(0, 1)

sns.lineplot(
    data=hyperparameters[hyperparameters['parameter'] == 'temperature'],
    x='iteration',
    y='mean',
    hue='chain',
    marker='o',
    palette='tab10',
    ax=axes[1]
)
axes[1].set_title('Temperature Hyperparameter Mean by Iteration')
axes[1].set_xlabel('Iteration')
axes[1].set_ylabel('Mean (temperature)')
axes[1].legend(title='Chain')
axes[1].set_ylim(0, 10)


plt.tight_layout()

hyperparameters[hyperparameters["lme"] == hyperparameters["lme"].max()]

from cpm.optimisation import minimise, FminBound


# Set up the fitting procedure
fit = FminBound(
    model=generative_model,  # Wrapper class with the model we specified from before
    data=data.groupby('ppt'),  # the data as a list of dictionaries
    minimisation=minimise.LogLikelihood.bernoulli,
    parallel=True,
    libraries=["numpy", "cpm", "pandas"],
    prior=True,
    ppt_identifier="ppt",
    display=False,
    number_of_starts=5,
    # everything below is optional and passed directly to the scipy implementation of the optimiser
    approx_grad=True

)

fit.optimise()
parameters_hierarchical = fit.export()
parameters_hierarchical.rename(
    columns={
        "x_0": "alpha",
        "x_1": "temperature",
    },
    inplace=True
)
parameters_hierarchical.head()

# Set up the fitting procedure
fit_no_prior = FminBound(
    model=generative_model,  # Wrapper class with the model we specified from before
    data=data.groupby('ppt'),  # the data as a list of dictionaries
    minimisation=minimise.LogLikelihood.bernoulli,
    parallel=True,
    libraries=["numpy", "cpm", "pandas"],
    prior=False,  # no prior
    ppt_identifier="ppt",
    display=False,
    number_of_starts=5,
    # everything below is optional and passed directly to the scipy implementation of the optimiser
    approx_grad=True

)

fit_no_prior.optimise()
parameters_no_prior = fit_no_prior.export()
parameters_no_prior.rename(
    columns={
        "x_0": "alpha",
        "x_1": "temperature",
    },
    inplace=True
)
parameters_no_prior.head()

from scipy.stats import truncnorm

best_hyperparameters = hyperparameters[hyperparameters["lme"] == hyperparameters["lme"].max()]
alpha_mean = best_hyperparameters[best_hyperparameters["parameter"] == "alpha"]["mean"].values[0]
alpha_sd = best_hyperparameters[best_hyperparameters["parameter"] == "alpha"]["sd"].values[0]
temperature_mean = best_hyperparameters[best_hyperparameters["parameter"] == "temperature"]["mean"].values[0]
temperature_sd = best_hyperparameters[best_hyperparameters["parameter"] == "temperature"]["sd"].values[0]

num_samples = 5000

fig =  plt.figure(figsize=(10, 10))

numpy.random.seed(984777324)
# For alpha: lower=0, upper=1
a_alpha = (0 - alpha_mean) / alpha_sd
b_alpha = (1 - alpha_mean) / alpha_sd
samples_alpha = truncnorm.rvs(a=a_alpha, b=b_alpha, loc=alpha_mean, scale=alpha_sd, size=num_samples)

# For temperature: lower=0, upper=10
a_temp = (0 - temperature_mean) / temperature_sd
b_temp = (10 - temperature_mean) / temperature_sd
samples_temperature = truncnorm.rvs(a=a_temp, b=b_temp, loc=temperature_mean, scale=temperature_sd, size=num_samples)

plt.scatter(x=parameters_hierarchical.alpha, y=parameters_hierarchical.temperature, color='red', label='Fitted Parameters', s=150, zorder=4, marker='+')
plt.scatter(x=parameters_no_prior.alpha, y=parameters_no_prior.temperature, color='black', label='Fitted Parameters', s=100, zorder=3, marker='o', alpha=0.5)
sns.kdeplot(x=samples_alpha, y=samples_temperature, cmap="Blues", thresh=0.05, linewidths=3)

plt.xlabel('Alpha')
plt.ylabel('Temperature')
plt.title('Comparison of Fitted Parameters and Hierarchical Priors')
plt.legend(['Hierarchical fit', 'Non-hierarchical fit', 'Hierarchical prior (KDE)'])
plt.tight_layout()