<jupyter_start><jupyter_text>Goals: 1. Find maximum-likelihood parameters for a model using gradient descent algorithms (based on SSE or LL)2. Model comparison using AIC and BIC3. Model and parameter recovery  Load packages<jupyter_code>import os, sys, glob, scipy
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
base_dir = os.path.realpath('.')
print(base_dir)<jupyter_output>/Users/jvanbaar/Dropbox (Brown)/Postdoc FHL/Teaching/Computational_modeling_tutorial<jupyter_text>Writing choice rules<jupyter_code>def greedy(utilities, params = []):
    
    choice = np.where(utilities == np.max(utilities))[0][0]
    
    return choice
def e_greedy(utilities, epsilon):
    
    random_number = np.random.rand()
    
    if random_number < epsilon:
        # Choose randomly
        choice = np.random.randint(0,2)
    else:
        # Choose deterministically
#         choice = np.where(utilites == np.max(utilities))
        choice = greedy(utilities)
    
    return choice
def softmax_with_choice(utilities, beta):
    # Beta is sometimes called 'inverse temperature' or 'explore parameter'
    
    # Version for choice A only => then prob_B will be 1 - prob_A
    prob_A = (np.exp(utilities[0] * beta) /
              (np.exp(utilities[0] * beta) + np.exp(utilities[1] * beta)))
    
    # General version (this works too, try it out!)
    probs = (np.exp(np.multiply(utilities, beta)) /
             (np.sum(np.exp(np.multiply(utilities, beta)))))
    
    random_number = np.random.rand()
    if random_number < prob_A:
        choice = 0
    else:
        choice = 1
    
    return choice
def softmax(utilities, beta):

    probs = (np.exp(np.multiply(utilities, beta)) /
             (np.sum(np.exp(np.multiply(utilities, beta)))))
    
    return probs<jupyter_output><empty_output><jupyter_text>Define our inequity aversion model again:<jupyter_code>def inequity_aversion(offer, alpha):
    
    # Inequity aversion model with one inequity term (= symmetric) and a money term
    payoffs = np.array([offer, 0])
    inequities = np.array([(100-offer)-offer, 0])
    utilities = payoffs - np.multiply(alpha,inequities)

    return utilities<jupyter_output><empty_output><jupyter_text>Simulate data for inequity aversion + softmax with varying parameter values<jupyter_code>def simulate_data(sub_data, alpha, beta, evaluate_choices = False):
    simulated_data = pd.DataFrame()
    for trial in sub_data['trial'].unique():
        trial_data = sub_data.query('trial == @trial')
        offer = trial_data['offer'].values[0]
        
        # Value function: inequity aversion ==> utilities
        utilities = inequity_aversion(offer, alpha)
        
        # Choice rule: softmax ==> probabilities
        choice_probabilities = softmax(utilities, beta)
        
        # Make a choice
        if evaluate_choices:
            choice = 0 if np.random.rand() < choice_probabilities[0] else 1
        
        # Store results
        if evaluate_choices:
            simulated_data = simulated_data.append(pd.DataFrame(
                [[trial, offer, choice_probabilities[0], choice]],
                columns=['trial', 'offer', 'probability_accept', 'choice'])).reset_index(drop=True)
        else:
            simulated_data = simulated_data.append(pd.DataFrame(
                [[trial, offer, choice_probabilities[0]]],
                columns=['trial', 'offer', 'probability_accept'])).reset_index(drop=True)
    return simulated_data<jupyter_output><empty_output><jupyter_text>Using the SciPy package to fit the model for us 1. Load Joey's data again2. Tell the SciPy.Optimise. to find the best fitting parameters  Load some data<jupyter_code>data = pd.read_csv(base_dir + '/Data/Study1_UG.csv')
data = data[['sub','trial','unfairness','choice']]
data['offer'] = 100 - data['unfairness']
data.head()<jupyter_output><empty_output><jupyter_text>Single subject:<jupyter_code>sub = 2
sub_data = data.query('sub == @sub').reset_index(drop=True)
sub_data.head()
plt.scatter(sub_data['unfairness'], sub_data['choice'])<jupyter_output><empty_output><jupyter_text>Implement the optimization routine<jupyter_code>def compute_negloglike(params, data):
    
    alpha = params[0]
    beta = params[1]
    
    # For each trial, compute the probability of observing THIS PARTICULAR BEHAVIOR given alpha and beta
    trial_probs = []
    for trial in data['trial'].unique():
        trial_data = data.query('trial == @trial')
        offer = trial_data['offer'].values[0]
        choice = trial_data['choice'].values[0]
        utilities = inequity_aversion(offer, alpha)
        choice_probabilities = softmax(utilities, beta)
        prob_of_observed = choice_probabilities[choice]
        trial_probs.append(prob_of_observed)
    
    # Multiply the probabilities of all trials to get the likelihood of the whole dataset –– here use the log trick!
    loglike = np.sum(np.log(trial_probs))
    
    return -loglike
import tqdm
n_iter = 10
objective_function = lambda params: compute_negloglike(params, sub_data)
fit_results = pd.DataFrame()
for i in tqdm.tqdm(range(n_iter)):
    x0 = np.random.rand(2)
    out = scipy.optimize.fmin(objective_function, x0, full_output=True)
    best_params, best_NLL = out[:2]
    fit_results = fit_results.append(pd.DataFrame([[i, best_params, best_NLL]],
                                     columns = ['iteration','best_params','best_NLL']))
fit_results.sort_values(by='best_NLL')<jupyter_output><empty_output><jupyter_text>Say you want to send this data to the GitHub repository for others to see<jupyter_code>fit_results.to_csv(base_dir + '/Results/fit_onesub_jeroen.csv')<jupyter_output><empty_output><jupyter_text>Now use git command line or GitHub Desktop to push this file to the repo (outside of notebook)<jupyter_code>!ls Results<jupyter_output>Study1_UG_Fits_Joey.csv fit_onesub_jeroen.csv<jupyter_text>Find the best-fitting parameters for all subjects using scipy.optimize.fmin Homework:1. Fit model for all subjects - find a good number of iterations to try2. Plot the best-fitting parameters for all subjects3. Do parameter recovery: a) pick 10 sets of parameters (alpha & beta), b) simulate data for each set, c) fit the model to the simulated data, d) repeat b and c 100 times, e) compare the 'recovered' parameters to the ones you picked in step a).  Load Joey's model fit results<jupyter_code>model_fits_allsubs = pd.read_csv(base_dir + '/Results/Study1_UG_Fits_Joey.csv',index_col=0)
model_fits_allsubs.head()<jupyter_output><empty_output><jupyter_text>Visualize the fitted parameters:<jupyter_code>sns.scatterplot(x = model_fits_allsubs.query('alpha < 10')['alpha'],
               y = model_fits_allsubs.query('alpha < 10')['beta'])
sns.distplot(model_fits_allsubs['alpha'], bins = 50)
sns.distplot(model_fits_allsubs['beta'], bins = 50)
sns.distplot(model_fits_allsubs.query('log_lik < -.1')['log_lik'])<jupyter_output><empty_output><jupyter_text>Compare to Lotte's results<jupyter_code>model_fits_lotte = pd.read_csv(base_dir + '/Results/Model_fit_results_Lotte.csv',
                              index_col = 0)
model_fits_lotte = model_fits_lotte.rename(columns = {'subject':'sub','tot_best_NLL':'NLL'})
model_fits_lotte['log_lik'] = -model_fits_lotte['NLL']
model_fits_lotte.head()
model_fits_joeylotte = model_fits_allsubs.merge(model_fits_lotte,
                           on = 'sub', suffixes = ['_joey','_lotte'])
model_fits_joeylotte.head()
sns.scatterplot(model_fits_joeylotte['alpha_joey'],
                model_fits_joeylotte['alpha_lotte'])
sns.scatterplot(model_fits_joeylotte['beta_joey'],
                model_fits_joeylotte['beta_lotte'])<jupyter_output><empty_output><jupyter_text>This looks pretty good: it seems like the subjects with alpha > 0 gave consistent results, while subjects with alpha < 0 gave erratic estimates for alpha. This means that alpha is not as well interpretable with < 0 values, which we can also conclude by reasoning about the model since all different alpha < 0 will yield the same result where the subject always accepts an UG offer. Actually going back to the original Fehr & Schmidt (1999) paper, where the IA model is proposed, they say alpha must be >= 0. So if we wanted, we could run the model again with that constraint.  Summary: how to deal with unreliable parameter fits? 1. Run the model more times to avoid local maxima.2. Apply boundaries to the model parameters, based on the theory behind the model or simulations of the model. E.g. Capping beta at 10, or writing the model such that alpha is between 0 and 1 with 0 meaning only caring about payoff and 1 meaning only caring about inequity.3. Run parameter recovery to test if the fits are reliable or not.  Parameter recovery  Logic: 1. Create a 'ground truth' model for which you KNOW the parameter values.2. Simulate data from that ground truth model.3. Fit your model to the simulated data to see if you can RECOVER the known parameters.<jupyter_code># Step 1: pick some ground truth alphas and betas
alpha = 1
beta = 2
# Step 2: simulate data
sim_data = simulate_data(sub_data, alpha, beta, evaluate_choices = True)
sim_data.head()
sns.scatterplot(data = sim_data, x = 'offer', y = 'choice')
# Step 3: fit model to these simulated data
n_iter = 10
objective_function = lambda params: compute_negloglike(params, sim_data)
fit_results = pd.DataFrame()
for i in range(n_iter):
    x0 = np.random.rand(2)
    out = scipy.optimize.fmin(objective_function, x0, full_output=True)
    best_params, best_NLL = out[:2]
    fit_results = fit_results.append(pd.DataFrame([[i, best_params[0], best_params[1], best_NLL]],
                                     columns = ['iteration','alpha','beta','best_NLL']))
best_NLL = np.min(fit_results['best_NLL'])
print(best_NLL)
best_solutions = fit_results.query('best_NLL == @best_NLL').copy()
best_solutions.head()
sns.scatterplot(data = best_solutions, x = 'alpha', y = 'beta')
plt.xlim([0,10]); plt.ylim([0,10]);
plt.scatter(alpha,beta)
best_solution = best_solutions.iloc[0]
best_solution
sim_best = simulate_data(sub_data, best_solution['alpha'], best_solution['beta'], evaluate_choices=True)
fig,ax = plt.subplots(nrows = 1, ncols = 2, figsize = [10,4])
ax[0].scatter(sim_data['offer'], sim_data['choice'])
ax[1].scatter(sim_best['offer'], sim_best['choice'])<jupyter_output><empty_output><jupyter_text>The parameters don't fit exactly but the data is the same... we therefore need more data in our experiment to get a more precise match of recovered to true parameters!<jupyter_code># Step 2: simulate data
sim_data = pd.DataFrame()
for n_experiments in range(100):
    to_append = simulate_data(sub_data, alpha, beta, evaluate_choices = True)
    sim_data = sim_data.append(to_append)
sim_data.head()
# Step 3: fit model to these simulated data
n_iter = 10
objective_function = lambda params: compute_negloglike(params, sim_data)
fit_results = pd.DataFrame()
for i in range(n_iter):
    x0 = np.random.rand(2)
    out = scipy.optimize.fmin(objective_function, x0, full_output=True)
    best_params, best_NLL = out[:2]
    fit_results = fit_results.append(pd.DataFrame([[i, best_params[0], best_params[1], best_NLL]],
                                     columns = ['iteration','alpha','beta','best_NLL']))
best_solutions = fit_results.query('best_NLL == @best_NLL').copy()
sns.scatterplot(data = best_solutions, x = 'alpha', y = 'beta')
plt.xlim([0,10]); plt.ylim([0,10]);
plt.scatter(alpha,beta)
sim_data_agg = sim_data.groupby(['offer'], as_index=False).mean()
sns.scatterplot(data = sim_data_agg, x = 'offer', y = 'choice')<jupyter_output><empty_output><jupyter_text>Let's try other alphas and betas<jupyter_code>alphas = np.arange(0,2,.4)
betas = np.arange(0,1,.2)
alphas, betas
n_iter = 1
recovery_results = pd.DataFrame()
for alpha in alphas:
    print(alpha)
    for beta in betas:
        sim_data = simulate_data(sub_data, alpha, beta, evaluate_choices=True)
        # Fit the model
        objective_function = lambda params: compute_negloglike(params, sim_data)
        fit_results = pd.DataFrame()
        for i in range(n_iter):
            x0 = np.random.rand(2)
            out = scipy.optimize.fmin(objective_function, x0, full_output=True)
            best_params, best_NLL = out[:2]
            fit_results = fit_results.append(pd.DataFrame([[i, best_params[0], best_params[1], best_NLL]],
                                             columns = ['iteration','alpha','beta','best_NLL']))
        to_append = pd.DataFrame([[alpha, beta, best_params[0], best_params[1], best_NLL]],
                                 columns = ['true_alpha','true_beta','recov_alpha','recov_beta','NLL'])
        recovery_results = recovery_results.append(to_append)
recovery_results.head()
sns.scatterplot(data = recovery_results.query('recov_alpha < 2 and true_alpha > 0'),
                x = 'true_alpha', y = 'recov_alpha')
plt.plot([0,1],[0,1],'k:')
sns.scatterplot(data = recovery_results.query('recov_alpha < 2 and true_alpha > 0'),
                x = 'true_beta', y = 'recov_beta')
plt.plot([0,10],[0,10],'k:')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11450148752765275, 'fraction_numerical': 0.012052788160805553, 'mean_word_length': 4.123094958968347, 'pattern_counts': {'":': 0, '<': 70, '<?xml version=': 0, '>': 66, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '526187', 'license': 'MIT', 'n_tokens_mistral': 4437, 'n_tokens_neox': 4215, 'n_words': 1071, 'path': '4.Automating_model_fitting_and_model_comparison.ipynb', 'repo_name': 'psychNerdJae/Computational_modeling_tutorial', 'size': '13109'}