<jupyter_start><jupyter_text>Psytrack Tutorial Link to paper: http://pillowlab.princeton.edu/pubs/Roy18_NeurIPS_dynamicPsychophys.pdfLast updated: February 16, 2020Psytrack version: 1.2.0---  Initialization<jupyter_code>%matplotlib inline
%load_ext autoreload
%autoreload 2

import numpy as np
from matplotlib import pyplot as plt
plt.rcParams['figure.dpi'] = 140

import psytrack as psy<jupyter_output><empty_output><jupyter_text>Quick DemonstrationThis is a quick, 2 minute _demonstration_ of Psytrack with a simulated dataset.\*\* **Please find a much more in-depth tutorial using real data below** \*\*--- Generate simulated dataThis includes generating psychometric weights ${\bf w}$, as well as the inputs ${\bf x}$ (sampled from a standard normal) and choices ${\bf y}$ (sampled according to the model using the weights and inputs).<jupyter_code>seed = 31
num_weights = 4
num_trials = 5000
hyper = {'sigma'   : 2**np.array([-4.0,-5.0,-6.0,-7.0]),
         'sigInit' : 2**np.array([ 0.0, 0.0, 0.0, 0.0])}

# Simulate
simData = psy.generateSim(K=num_weights, N=num_trials, hyper=hyper,
                          boundary=6.0, iterations=1, seed=seed, savePath=None)

# Plot
psy.plot_weights(simData['W'].T);
plt.ylim(-3.6,3.6);<jupyter_output><empty_output><jupyter_text>Recover the weights from the simulated behaviorThat is, given inputs ${\bf x}$ and choices ${\bf y}$, recover the psychometric weights ${\bf w}$._Note: This takes approximately 60 seconds_<jupyter_code>rec = psy.recoverSim(simData)<jupyter_output><empty_output><jupyter_text>Plot the true weights from above (now in black) with the recovered weights (in color)Recovered weights also have shading to indicate a 95\% credible interval.<jupyter_code>psy.plot_weights(rec['wMode'], errorbar=rec["hess_info"]["W_std"])
plt.plot(simData['W'], c="k", ls="-", alpha=0.5, lw=0.75, zorder=0)
plt.ylim(-3.6,3.6);<jupyter_output><empty_output><jupyter_text>Plot recovered smoothness hyperparameters $\sigma_k$ over the true hyperparameters (black lines)Recovered hyperparameters plotted with $\pm$2SE bars.<jupyter_code>true_sigma = np.log2(rec['input']['sigma'])
avg_sigma = np.log2(rec['hyp']['sigma'])
err_sigma = rec['hess_info']['hyp_std']

plt.figure(figsize=(2,2))
colors = np.unique(list(psy.COLORS.values()))
for i in range(num_weights):
    plt.plot(i, true_sigma[i], color="black", marker="_", markersize=12, zorder=0)
    plt.errorbar([i], avg_sigma[i], yerr=2*err_sigma[i], color=colors[i], lw=1, marker='o', markersize=5)

plt.xticks([0,1,2,3]); plt.yticks(np.arange(-8,-2))
plt.gca().set_xticklabels([r"$\sigma_1$", r"$\sigma_2$", r"$\sigma_3$", r"$\sigma_4$"])
plt.xlim(-0.5,3.5); plt.ylim(-7.5,-3.5)
plt.gca().spines['right'].set_visible(False)
plt.gca().spines['top'].set_visible(False)
plt.ylabel(r"$\log_2(\sigma)$");<jupyter_output><empty_output><jupyter_text>------ Full Tutorial with Real DataDatasets handled by Psytrack are specific to an individual animal and are stored as a Python dictionary<jupyter_code># Extract premade dataset from npz
D = np.load('sampleRatData.npz', allow_pickle=True)['D'].item()

print("The keys of the dict for this example animal:\n   ", list(D.keys()))<jupyter_output><empty_output><jupyter_text>Of these keys, only `y` and `inputs` are **required** for analysis of the dataset, all others are *optional*.---`y` should be a 1D array of the animal's choice on each trial. Currently, the analysis only works for two-alternative forced choice tasks, and so there should only be two options on each trial (error or omission trials are typically discarded from the analysis).The two options (A or B, Left or Right, etc.)  must be mapped to {1, 2} _or_ {0, 1} in `y` (not -1 and +1). In this example, Left=1 and Right=2.<jupyter_code>print("The shape of y:   ", D['y'].shape)
print("The number of trials:   N =", D['y'].shape[0])
print("The unique entries of y:   ", np.unique(D['y']))<jupyter_output><empty_output><jupyter_text>---`inputs` is itself another dictionary, containing arbitrary keys. Each of these keys represents a _potential_ input into the model and must be a 2D array of shape $(N, M)$ where $N$ is the number of trials. The number of columns $M$ is arbitrary, and the $i^{th}$ column is typically used to encode information from $i$ time steps previous.For example, in our example data set the key `s1` encodes the (normalized) stimulus values heard on each trial. `s1[7,0]` would encode the stimulus heard on the 7th trial where as both `s1[6,0]` and `s1[7,1]` would encode the stimulus heard on the 6th trial. The information is redundant, but allows for all feasible regressors to predicting behavior of trial $i$ to be accessible by referencing the $i^{th}$ row of the respective input array.<jupyter_code>print("The keys of inputs:\n   ", list(D['inputs'].keys()))

print("\nThe shape of s1:", D['inputs']['s1'].shape)
print("s1[7]   : ", D['inputs']['s1'][7])
print("s1[6,0] : ", D['inputs']['s1'][6,0])
print("s1[7,1] : ", D['inputs']['s1'][7,1])<jupyter_output><empty_output><jupyter_text>---Other keys are for convenience and are _optional_ : `name` stores the name of the animal, `answer` is an easy reference as to what the correct choice was on a given trial, and `correct` is an easy reference as to whether the animal made the correct choice on a given trial. The model only needs to know what the animal _actually_ did, not what the animal _ought_ to have done!`dayLength` is an array storing the number of trials that occurred in each session of training. Taking a cumulative sum will give you the indices at which each new session begins. This is **not** optional for the analysis if one wishes to use the `sigmaDay` functionality (see Section 3.3 in paper). --- Fitting the dataOnce you have your data arranged in the proper format, you can now run the analysis! The fitting function is called `hyperOpt()` and before using it, you must decide on 3 inputs:   1) `weights` : which of your inputs should you fit.   2) `hyper` : what hyperparameters should your model have and how should they be initialized.   3) `optList` : what subset of the hyperparameters should be optimized. `weights` is a dictionary where the keys correspond to the keys in your dataset's `inputs` dictionary; the key values are a non-negative integer indicating how many of the columns of that value in `inputs` should be used for fitting the model, where each included column corresponds to a new weight. You can also include in `weights` the special key `bias` which need not be included in `inputs` --- this will simply create an input of all 1s.<jupyter_code>weights = {'bias': 1,  # a special key
           's1': 1,    # use only the first column of s1 from inputs
           's2': 1}    # use only the first column of s2 from inputs

# It is often useful to have the total number of weights K in your model
K = np.sum([weights[i] for i in weights.keys()])<jupyter_output><empty_output><jupyter_text>`hyper` is a dictionary that indicates what hyperparameters your model will have. There are 3 types:1) `sigma` : required, controls trial-to-trial variability.2) `sigInit` : optional, controls the variability on the very first trial (e.g. how close weights must initialize to 0). It is often best to include this hyperparameter and set it to a high value, as you often prefer your data to determine where the weights ought to initialize. Otherwise, `sigInit` will be set equal to `sigma`.3) `sigDay` : optional, controls variability between sessions (e.g. between the last trial of one session and the first trial of the next session). If this key is included, then your dataset must also have the key `dayLength` (as described above), to indicate the trials where the `sigDay` variability should supercede the standard `sigma` variability.For each hyperparameter key included in the `hyper` dictionary, the corresponding value is the initial value of the hyperparameter. If you are optimizing over a particular hyperparameter (see the `optList` setting below), than the initial value is not so important as the fitting procedure will eventually converge to the optimal setting. However, if you are *not* optimizing, then the initial value set will be the *fixed* value of the hyperparameter.Finally, for each hyperparameter key in `hyper`, you must specify your initializations as a 1D array with length $K$. If you instead provide only a single value, then the optimizer will assume that you want the same hyperparameter to apply to every weight (as opposed to each weight having it's own).<jupyter_code>hyper= {'sigInit': 2**4.,      # Set to a single, large value for all weights. Will not be optimized further.
        'sigma': [2**-4.]*K,   # Each weight will have it's own sigma optimized, but all are initialized the same
        'sigDay': None}        # Indicates that session boundaries will be ignored in the optimization<jupyter_output><empty_output><jupyter_text>`optList` is a list of the subset of hyperparameters in `hyper` that you wish to optimize over in your model. It is typically unnecessary to optimize over `sigInit` -- a single, large, fixed value is usually best. Thus, there are 4 reasonable ways to specify `optList`:1) `optList = []` : this means that none of the hyperparameters in `hyper` will be optimized. The optimizer will find the best set of weight trajectories given the fixed hyperparameter values specified.2) `optList = ['sigma']` : only the `sigma` hyperparameter will be optimized (as we will do below). If `sigDay` is included in `hyper`, then this means that the model _will_ include `sigDay`, but that the initial value specified will not be optimized.3) `optList = ['sigDay']` : only the `sigDay` hyperparameter will be optimized, while the `sigma` hyperparameter remains fixed.4) `optList = ['sigma', 'sigDay']` : both the `sigma` and `sigDay` hyperparameters will be optimized.<jupyter_code>optList = ['sigma']<jupyter_output><empty_output><jupyter_text>Now that we have specified `weights`, `hyper`, and `optList`, we can fit our model with `hyperOpt()`! The function takes your dataset `D` plus the three additional inputs, and returns 4 things:1) `hyp` : a dictionary of the optimized hyperparameters2) `evd` : the approximate log-evidence of the optimized model3) `wMode` : the weight trajectories of the optimized model4) `hess_info` : a dictionary of sparse terms that relate to the Hessian of the optimal model. By default, this also includes the posterior credible intervals on the weights, under the key `W_std`. This behavior can be altered by changing the optional argument `hess_calc` in `hyperOpt()` (see function documentation for more details).Run times will depend on the number of trials $N$ and weights $K$, as well as the number of hyperparameters being fit. To speed things up a bit, we will use the `trim` function to shrink our dataset of 20K trials to just the first 10K._Note: this should take < 60 seconds._<jupyter_code>new_D = psy.trim(D, END=10000)  # trim dataset to first 10,000 trials
hyp, evd, wMode, hess_info = psy.hyperOpt(new_D, hyper, weights, optList)<jupyter_output><empty_output><jupyter_text>--- Visualizing the resultsPsytrack includes a few plotting functions for visualizing the results of the fit. To see the weight trajectories, use `plot_weights`.<jupyter_code>fig = psy.plot_weights(wMode, weights)<jupyter_output><empty_output><jupyter_text>Adding credible intervals on the weights can be done with the `errorbar` keyword argument, and vertical lines indicating session boundaries can be added with `days`. Adjustments to the resulting plot can be made by editing the figure returned by the function directly.<jupyter_code>fig = psy.plot_weights(wMode, weights, days=new_D["dayLength"], errorbar=hess_info["W_std"])<jupyter_output><empty_output><jupyter_text>We can also generate two additional plots with useful information: 1) A performance plot, tracking the animal's task accuracy, smoothed over trials2) A bias plot, tracking the animal's choice bias, smoothed over trialsBoth of these plots will calculate their respective values directly from the data, with 2SD error bars.<jupyter_code>fig_perf = psy.plot_performance(new_D)
fig_bias = psy.plot_bias(new_D)<jupyter_output><empty_output><jupyter_text>--- Cross-validationThe performance and bias plots above use the data directly to plot the *empirical* acuracy and bias of the animal. We can use our model to make predictions about the animal's accuracy and bias using the weight trajectories found by the model, to see if there is agreement. However, to make true predictions, we need to make predictions on trials that were held-out from the fitting procedure.We can do this using the built-in cross-validation functions `crossValidate`. This is also useful if you'd like to compare different models via cross-validated log-likelihood, rather than approximate model evidence.Similar to `hyperOpt`, `crossValidate` receives the same inputs as well an additional input `F`, which controls the number of folds in the cross-validation procedure (note that `F` must cleanly divide the number of trials in the dataset `N` -- the `trim` function can be used to make this work). `crossValidate` will then divide the dataset into `F` training and testing datasets, fitting each of the `F` training data sets. It will return `xval_logli`, the total cross-validated log-likelihood of the `F` test sets, as well `xval_pL`, the cross-validated $P(y=0)$ for each trial._Note: Since we're fitting `F` models, this can be fairly time consuming. This particular example should take about 10 minutes._<jupyter_code>xval_logli, xval_pL = psy.crossValidate(new_D, hyper, weights, optList, F=10, seed=41)<jupyter_output><empty_output><jupyter_text>Using the cross-validated prediction that the animal will go Left on each trial, `xval_pL`, we can overlay the plots of the empirical performance and bias above with a line showing the model's prediction. Fortunately, we see a close agreement!<jupyter_code>fig_perf_xval = psy.plot_performance(new_D, xval_pL=xval_pL)
fig_bias_xval = psy.plot_bias(new_D, xval_pL=xval_pL)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09832813828279965, 'fraction_numerical': 0.013672428449985832, 'mean_word_length': 5.257535460992908, 'pattern_counts': {'":': 0, '<': 70, '<?xml version=': 0, '>': 69, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '247539', 'license': 'MIT', 'n_tokens_mistral': 4177, 'n_tokens_neox': 3984, 'n_words': 1800, 'path': 'psytrack/examples/ExampleNotebook.ipynb', 'repo_name': 'andialbert/psytrack', 'size': '14116'}