<jupyter_start><jupyter_text>Model comparison<NAME>, March 2022<jupyter_code>%matplotlib inline

import os, sys
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle

import seaborn as sns
import teaching_models as teach
from os.path import join as opj
from scipy.spatial import distance
from ast import literal_eval

sys.path.append('..')
from utils import read_json, gsearch, str_extract, int_extract, print_list

sns.set_style('white')
sns.set_context('talk')<jupyter_output>Loading teaching problems
{'A': [[0, 0, 1, 1, 0, 0],
       [0, 1, 1, 1, 1, 0],
       [1, 1, 1, 1, 1, 1],
       [1, 1, 1, 1, 1, 1],
       [0, 1, 1, 1, 1, 0],
       [0, 0, 1, 1, 0, 0]],
 'B': [[1, 1, 1, 0, 0, 0],
       [1, 1, 1, 0, 0, 0],
       [1, 1, 1, 1, 0, 0],
       [0, 0, 1, 1, 0, 0],
       [0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0]],
 'C': [[0, 0, 0, 1, 1, 1],
       [0, 0, 0, 1, 1, 1],
       [0, 0, 1, 1, 1, 1],
       [0, 0, 1, 1, 0, 0],
       [0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0]],
 'D': [[0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0],
       [0, 0, 1, 1, 0, 0],
       [1, 1, 1, 1, 1, 1],
       [1, 1, 1, 1, 1, 1],
       [1, 1, 1, 1, 1, 1]]}
Participants excluded from analysis:
[3, 17]
Cleaning up human behavioral data
    subject  run  block_idx  ex_idx first_movement   start  problem  \
0         1    1          0       0         [0, 1]  [0, 0]       22   
1         1    1          0       1         [0, 1]  [0, 0]       22   
3         1    1          1   [...]<jupyter_text>Model recovery Find model recovery files:<jupyter_code>model_files = gsearch('outputs/simulated_model_comparison/*pxp.txt')
model_files.sort()
print_list(model_files)<jupyter_output>Found 11 items
outputs/simulated_model_comparison/model-cost_pxp.txt
outputs/simulated_model_comparison/model-literal_cost_pxp.txt
outputs/simulated_model_comparison/model-literal_pref_cost_pxp.txt
outputs/simulated_model_comparison/model-literal_pref_pxp.txt
outputs/simulated_model_comparison/model-literal_pxp.txt
outputs/simulated_model_comparison/model-pragmatic_cost_pxp.txt
outputs/simulated_model_comparison/model-pragmatic_pref_cost_pxp.txt
outputs/simulated_model_comparison/model-pragmatic_pref_pxp.txt
outputs/simulated_model_comparison/model-pragmatic_pxp.txt
outputs/simulated_model_comparison/model-pref_cost_pxp.txt
...<jupyter_text>Data wrangling<jupyter_code>pxp_list = []
model_levels = ['pragmatic_pref_cost',  'pragmatic_pref',  'pragmatic_cost', 'pragmatic', 'literal_pref_cost', 'literal_pref',  'literal_cost', 'literal', 'pref_cost',  'pref', 'cost']
for f in model_files:
    pxp = np.loadtxt(f,delimiter=',')

    label_f = f.replace('pxp', 'model_labels')
    labels = np.loadtxt(label_f, dtype=str)

    f_df = pd.DataFrame(list(zip(labels, pxp)), columns=['fit_model', 'pxp'])
    f_df['true_model'] = str_extract('(?<=model-).+(?=_pxp)', f)
    
    pxp_list.append(f_df)
    
pxp_df = pd.concat(pxp_list)
pxp_df['fit_model'] = pxp_df.fit_model.astype('category').cat.reorder_categories(model_levels)
pxp_df['true_model'] = pxp_df.true_model.astype('category').cat.reorder_categories(model_levels)

print(pxp_df.shape)
pxp_df.head()<jupyter_output>(121, 3)<jupyter_text>Make confusion matrix<jupyter_code>confusion_mtx = pxp_df.pivot(index='fit_model', columns='true_model', values='pxp')
ax = sns.heatmap(confusion_mtx, cmap='viridis', square=True)
ax.set(xlabel='True model', ylabel='Predicted model')
plt.savefig('plots/model_recovery.png', bbox_inches='tight')<jupyter_output><empty_output><jupyter_text>Model comparison: Which model best captures human behavior? Read model evidence file:<jupyter_code>human_pxp = np.loadtxt('outputs/human_model_comparison/human_model_pxp.txt', delimiter=',')
human_labels = np.loadtxt('outputs/human_model_comparison/model_labels.txt', dtype=str)

human_df = pd.DataFrame({'model': human_labels, 'pxp': human_pxp})
human_df['model'] = human_df['model'].astype('category').cat.reorder_categories(model_levels + ['strong'])
human_df = human_df.sort_values(by='model').reset_index(drop=True)

human_df<jupyter_output><empty_output><jupyter_text>Read probability that each model generated each teacher:<jupyter_code>fig,ax=plt.subplots(figsize=(12,4.5))
sns.barplot(data=human_df, x='model', y='pxp', ax=ax, fc='#333')
ax.set(xlabel='', ylabel='Protected exceedance\nprobabilities')
ax.set_xticklabels(ax.get_xticklabels(),rotation = 90)
plt.savefig('plots/model_comparison_results.pdf', bbox_inches='tight')<jupyter_output><empty_output><jupyter_text>Predictive checks: How well does the winning model capture qualitative patterns in human behavior? Helper function: Return frequency of responses as a 6x6 matrix for plotting<jupyter_code>def example_matrix(group):
    freq = np.zeros(36)
    freq[group.example] = group.freq
    freq_mtx = np.reshape(freq, (6,6))
    
    return freq_mtx<jupyter_output><empty_output><jupyter_text>Find simulation files<jupyter_code>sim_files = gsearch('outputs', 'simulated_from_human', '*pref_cost.csv')
sim_files.sort()

print_list(sim_files)<jupyter_output>Found 3 items
outputs/simulated_from_human/simulated_data_model-literal_pref_cost.csv
outputs/simulated_from_human/simulated_data_model-pragmatic_pref_cost.csv
outputs/simulated_from_human/simulated_data_model-pref_cost.csv<jupyter_text>Wrangle simulation data<jupyter_code>sim_df = pd.concat([pd.read_csv(f) for f in sim_files])

print(sim_df.shape)
sim_df.head()<jupyter_output>(10080, 7)<jupyter_text>Plot average simulation results against human data<jupyter_code>def predictive_check(prob_idx):
    # Mask plots using true hypothesis
    prob = teach.problems[prob_idx]
    mask = 1-np.array(prob['A'])

    # Average human responses
    prob_human = (
        teach.human_df
        [teach.human_df.problem == prob_idx]
        .copy()
        .groupby('example')
        .agg({'subject': 'count'})
        .reset_index()
        .rename(columns={'subject': 'n'})
    )
    prob_human['freq'] = prob_human.n/prob_human.n.sum()

    # Average simulated responses
    prob_sim = (
        sim_df[sim_df.problem == prob_idx]
        .copy()
        .groupby(['model', 'example'])
        .agg({'subject': 'count'})
        .reset_index()
        .rename(columns={'subject': 'n'})
    )
    prob_sim['n_total'] = prob_sim.groupby(['model']).n.transform('sum')
    prob_sim['freq'] = prob_sim.n/prob_sim.n_total

    # Plot human data
    n_models = prob_sim.model.nunique()
    fig,axes = plt.subplots(1,n_models+1, figsize=(4*(n_models+1)*1.25, 4))
    human_mtx = example_matrix(prob_human)
    sns.heatmap(human_mtx, cmap='viridis', square=True, ax=axes[0], mask=mask, lw=2, linecolor='#fff')
    axes[0].set(title='Human', xticklabels=[], yticklabels=[])

    # Plot simulated data
    ax_counter = 1
    for name, group in prob_sim.groupby('model'):
        freq_mtx = example_matrix(group)
        sns.heatmap(freq_mtx, cmap='viridis', square=True, ax=axes[ax_counter], mask=mask, lw=1, linecolor='#fff')
        axes[ax_counter].set(title=name, xticklabels=[], yticklabels=[])

        ax_counter += 1
#_ = teach.plot_problem(3)
predictive_check(3)
plt.savefig('plots/predictive_check.png')<jupyter_output><empty_output><jupyter_text>And here's a counterexample where the pragmatic model does particularly well:<jupyter_code>_ = teach.plot_problem(11)
predictive_check(11)<jupyter_output><empty_output><jupyter_text>Sanity check: Plot idealized model predictions Informational value<jupyter_code>prob = teach.problems[3]
prob_mask = 1-np.array(prob['A'])

pD,pH = teach.strong_sampling(3)
info_value = np.log(pH)
info_mtx = teach.sampling_matrix(info_value)

ax = sns.heatmap(info_mtx, square=True, cbar=False, cmap='viridis', lw=1, linecolor='#fff', mask=prob_mask)
ax.set(xticklabels=[], yticklabels=[])

plt.savefig('plots/info_value.png')<jupyter_output><ipython-input-14-68d3f3c1fba0>:5: RuntimeWarning: divide by zero encountered in log
  info_value = np.log(pH)
/n/helmod/apps/centos7/Core/Anaconda3/2020.11-jupyterood-fasrc01/x/lib/python3.8/site-packages/matplotlib/colors.py:1062: RuntimeWarning: overflow encountered in true_divide
  resdat /= (vmax - vmin)<jupyter_text>Edge preference:<jupyter_code>prob = teach.problems[3]
prob_mask = 1-np.array(prob['A'])

edge_df = teach.edge_pref(prob)
edge_mtx = teach.sampling_matrix(edge_df)
ax = sns.heatmap(edge_mtx, square=True, cbar=False, cmap='viridis', lw=1, linecolor='#fff', mask=prob_mask)
ax.set(xticklabels=[], yticklabels=[])

plt.savefig('plots/edge_pref.png')<jupyter_output><empty_output><jupyter_text>Distance:<jupyter_code>cost_df = teach.movement_cost(3)
cost_mtx = teach.sampling_matrix(cost_df)

ax = sns.heatmap(cost_mtx, square=True, cbar=False, cmap='viridis', lw=1, linecolor='#fff', mask=prob_mask)
ax.set(xticklabels=[], yticklabels=[])
ax.add_patch(Rectangle((0, 0), 1.0, 1.0, edgecolor='red', fill=False, lw=2))

plt.savefig('plots/distance.png')<jupyter_output><empty_output><jupyter_text>Utility:<jupyter_code>pD, pH = teach.utility_sampling(3, sampling_fun=teach.strong_sampling, pref_fun=teach.edge_pref)
utility_mtx = teach.sampling_matrix(pD)

ax = sns.heatmap(utility_mtx, square=True, cbar=False, cmap='viridis', lw=1, linecolor='#fff', mask=prob_mask)
ax.set(xticklabels=[], yticklabels=[])
ax.add_patch(Rectangle((0, 0), 1.0, 1.0, edgecolor='red', fill=False, lw=2))

plt.savefig('plots/utility.png')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15051939792240832, 'fraction_numerical': 0.029255882976468094, 'mean_word_length': 5.1385816525699415, 'pattern_counts': {'":': 0, '<': 58, '<?xml version=': 0, '>': 57, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '77402', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3855, 'n_tokens_neox': 3578, 'n_words': 651, 'path': '2_behavioral/7_model_comparison.ipynb', 'repo_name': 'nataliavelez/teaching_fmri', 'size': '9441'}