Estimating Empirical Priors (cpm.hierarchical)¶

experiment['observed'] = experiment['responses'] - 1

from cpm.generators import Parameters, Value

## define our parameters
parameters = Parameters(
    # freely varying parameters are indicated by specifying priors
    alpha=Value(
        value=0.5,
        lower=1e-10,
        upper=1,
        prior="truncated_normal",
        args={"mean": 0.6, "sd": 0.20}, # we pick some random values here for the true prior
    ),
    temperature=Value(
        value=1,
        lower=0,
        upper=10,
        prior="truncated_normal",
        args={"mean": 2, "sd": 2.5}, # we pick some random values here for the true prior
    ),
    # everything without a prior is part of the initial state of the model or fixed constructs (e.g. exemplars in general-context models of categorization)
    values = numpy.array([0.25, 0.25, 0.25, 0.25]))

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

@ipp.require("numpy")
def model(parameters, trial):
    # pull out the parameters
    alpha = parameters.alpha
    temperature = parameters.temperature
    values = numpy.asarray(parameters.values)
    # pull out the trial information
    stimulus = numpy.array([trial.stimulus_left, trial.stimulus_right]).astype(int)
    feedback = numpy.array([trial.reward_left, trial.reward_right])
    human_choice = trial.responses.astype(int) - 1 # convert to 0-indexing
    mute = numpy.zeros(4)  # mute learning for all cues not presented

    # activate the value of each available action
    # here there are two possible actions, that can take up on 4 different values
    # so we subset the values to only include the ones that are activated...
    # ...according to which stimuli was presented
    activations = values[stimulus - 1]
    # calculate a policy based on the activations
    response = cpm.models.decision.Softmax(activations=activations, temperature=temperature)
    response.compute() # compute the policy
    if numpy.isnan(response.policies).any():
        # if the policy is NaN for a given action, then we need to set it to 1
        print(f"NaN in policy with parameters: {alpha.value}, {epsilon.value}\n")
        print(response.policies)
        response.policies[numpy.isnan(response.policies)] = 1

    model_choice = response.choice() # get the choice based on the policy
    reward = feedback[human_choice] # get the reward of the chosen action
    # update the value of the chosen action
    mute[stimulus[human_choice] - 1] = 1 # unmute the learning for the chosen action
    teacher = numpy.array([reward])
    update = cpm.models.learning.SeparableRule(weights=values, feedback=teacher, input=mute, alpha=alpha)
    update.compute()
    values += update.weights.flatten()
    ## compile output
    output = {
        "trial"    : trial.trial.astype(int),               # trial number
        "policy"   : response.policies,         # policies
        "response" : model_choice,              # choice based on the policy
        "reward"   : reward,                    # reward of the chosen action
        "values"   : values,                    # updated values
        "change"   : update.weights,            # change in the values
        "activation" : activations.flatten(),     # activation of the values
        "dependent"  : numpy.array([response.policies[1]]),        # dependent variable
    }
    return output

true = parameters.sample(experiment.ppt.unique().shape[0]) # sample a set of parameters

simulator = Simulator(wrapper=wrapper,
                      data=experiment.groupby("ppt"),
                      parameters=true,
                      ) # data grouped by participants
simulator.run()

generated = simulator.export()
experiment.response = generated.response

hyperparameters = eb.hyperparameters
hyper = hyperparameters[hyperparameters.lme == hyperparameters.lme.max()]
means = [hyper[hyper.parameter == "alpha"].squeeze()['mean'],
         hyper[hyper.parameter == "temperature"].squeeze()['mean']]

ssd = [hyper[hyper.parameter == "alpha"].squeeze()['sd'],
       hyper[hyper.parameter == "temperature"].squeeze()['sd']]

## redefine parameters
estimated_parameters = Parameters(
    # freely varying parameters are indicated by specifying priors
    alpha=Value(
        value=0.5,
        lower=1e-10,
        upper=1,
        prior="truncated_normal",
        args={"mean": means[0], "sd": ssd[0]},
    ),
    temperature=Value(
        value=1,
        lower=0,
        upper=10,
        prior="truncated_normal",
        args={"mean": means[1], "sd": ssd[1]},
    ),
    # everything without a prior is part of the initial state of the model or fixed constructs (e.g. exemplars in general-context models of categorization)
    values = numpy.array([0.25, 0.25, 0.25, 0.25]))

## reinitialise the model
wrapper_with_prior = Wrapper(model=model, parameters=estimated_parameters, data=experiment[experiment.ppt == 1])
fit_with_prior = FminBound(
    model=wrapper,  # Wrapper class with the model we specified from before
    data=experiment,  # the data as a list of dictionaries
    minimisation=minimise.LogLikelihood.bernoulli,
    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_with_prior.optimise()

original = pandas.DataFrame(true)
non_hierarchical = fit.export()
with_hierarchical = fit_with_prior.export()