Lab Assignments

Part 1: Matching Law

This week, we are focusing on the generalized matching equation (GME). The purpose of this lab is to use your programming skills to fit the GME to participant data and interpret the outputs. In the folder, there are 10 hypothetical datasets with participant data common to matching law experiments. Your job is what a researcher's job might be. That is, see how well the data are fit by the GME for each individual as well as any trends that might be worth discussing at the group level.

The following packages will likely help accomplish this task: pandas, numpy, matplotlib, seaborn, and scipy (linregress function). These are not the only ways you might accomplish the goal, but they certainly have everything you need.

Part 2: Discounting

We are also focusing on the family of equations within the discounting area of the literature. The purpose of this part of the lab is to use your programming skills to fit the hyperbolic and hyperboloid (Green & Myerson) equations and compute Area Under the Curve for participant data, and interpret the outputs.

In the folder, there is a single dataset with participant data common to discounting experiments. Specifically, each row has a set of indifference points relative to seven different indifference points (columns) along with information about the amount and commodity specific to those indifference points.

Your job is what a researcher's job might be. That is: (1) see how well the data are fit by the hyperbolic and hyperboloid models; (2) once fit, see if there are any trends in discounting relative to changes in amount, commodity, or sign (gain vs. loss).

The following packages will likely help accomplish this task: pandas, numpy, matplotlib, seaborn, scipy (for fitting and calculating AUC), and scikit-learn (for r2 values).

As a bonus, in the folder, you will also see a params_data.csv. These are the k and s parameters of the hyperboloid used to create each participant's data, before zero-mean Gaussian noise (SD 4% of the amount) was added at every delay except zero. You could compare how well your fits recover the original parameters when noise is involved.

Demand Lab

This week, we are focusing on the family of equations within the demand area of behavior analysis. The purpose of this lab is to use your programming skills to explore demand. In this assignment, you will fit the exponential demand equation proposed by Hursh & Silberberg (2008) to participant consumption data.

In this folder, there is a single dataset with participant data common in demand experiments. Specifically, each row has a single participant with consumption at a range of price points (i.e., $0.01,$0.05, $0.25,$1, $5,$10, $33,$100, $250,$500, $1000). The Q0 and alpha values used to generate each participant's data are also present. However, noise was added once generated. Your job is to calculate Q0 from the participant's data directly, calculate the k value from the dataset as a whole, and then systematically loop through each participant and fit the exponential equation, estimating alpha for each participant.

General Steps

1. Read in the dataset
2. Transform the dataset to be in long form
3. Transform the price to a numeric value
4. Calculate the parameter k (the range of nonzero consumption in log10 units, plus headroom) or choose its constant value
5. Create a function that calculates the consumption Q using the exponential demand equation in the chapter's form, with the exponent -alpha * Q0 * C
6. Apply the function to each row in the dataset to identify the predicted consumption
7. Graph and analyze the results.
8. Calculate the goodness of fit metrics and display them to interpret the model's goodness of fit

Things to Watch Out For

• The log of 0 is undefined because there's no number you can raise that will give you zero.
• Often to avoid the log(0) issue, people will add a constant to all values in the dataset before taking the log (i.e., all numbers in the series are adjusted up by 0.001 before the log transform is applied).
• Here, k is computed from the nonzero consumption values, and the model is fit in consumption units (predicted Q against observed Q), so zeros at the highest prices are ordinary observations rather than a problem.

Reference

Hursh, S. R., & Silberberg, A. (2008). Economic demand and essential value. Psychological Review, 115(1), 186-198. https://doi.org/10.1037/0033-295X.115.1.186

Associative Learning Models Lab

This week, we focus on respondent conditioning -- stimulus-stimulus learning -- and the quantitative models that describe how organisms learn predictive relationships between events. Both models in this lab are recursive: the state of the system at trial t (the associative strengths, and for Mackintosh the attention weights) becomes part of the input at trial t+1. Formally:

X(t+1) = f(X(t), I(t))

where X(t) is the state at time t, I(t) is the input or event at time t (a CS, the US, a prediction error), and f is the transition rule. Because both parts are simulations, no data files are required.

Part 1: The Rescorla-Wagner Model

The Rescorla-Wagner (1972) model is the foundational quantitative account of associative learning. Its central claim is that learning is driven by prediction error: associative strength changes only to the extent that the outcome of a trial differs from what the organism already expected. For each stimulus present on a trial,

ΔV_i = α_i * β * (λ - V_total), where V_total = Σ V_i

Here α_i is the salience of stimulus i, β is the learning rate of the US, λ is the asymptote supported by the US (1.0 when it occurs, 0.0 when it does not), and V_total is the summed associative strength of all stimuli present on that trial.

Your task is to implement this update rule as a function and use it to reproduce four classic phenomena that all emerge from the shared prediction-error term:

• Blocking: a CS that already predicts the US prevents learning to a newly added CS.
• Overshadowing: when two CSs are trained together, the more salient one captures more associative strength.
• Overexpectation: two separately conditioned CSs, then trained together, both lose strength because their summed prediction overshoots λ.
• Conditioned inhibition: a CS paired with the absence of an otherwise-predicted US acquires negative associative strength.

Hints:

• Represent each trial as the set of stimuli present plus the US asymptote λ for that trial. Build the phenomenon by choosing the right sequence of trials across phases.
• Compute V_total from the stimuli present before updating, then apply ΔV to each present stimulus.

Part 2: The Mackintosh Model (Recursive Associability)

The Mackintosh (1975) model extends prediction-error learning with a dynamic associability (attention) term. A stimulus that predicts the US better than its competitors gains associability across trials, while a poorer predictor loses it. This lets the model capture attentional phenomena that fixed-salience Rescorla-Wagner cannot.

Your task is to implement the Mackintosh update rules for association strength V and associability alpha, then write simulation functions for (a) basic conditioning, (b) overshadowing, and (c) blocking, and visualize how V and alpha evolve across trials for a range of learning rates (theta).

Hints:

• Recursive nature: both the associability (alpha) and the association strength (V) update based on their previous values.
• Phenomenon validation: your outputs should show basic conditioning (V rises toward λ while alpha adjusts), overshadowing (CS1 with higher alpha gains more V than CS2), and blocking (CS2 shows minimal learning when CS1 already predicts the US).

References

• Rescorla, R. A., & Wagner, A. R. (1972). A theory of Pavlovian conditioning: Variations in the effectiveness of reinforcement and nonreinforcement.
• Mackintosh, N. J. (1975). A theory of attention: Variations in the associability of stimuli with reinforcement. Psychological Review, 82(4), 276-298.

Behavioral Momentum and Response Persistence Lab

Behavioral momentum theory (Nevin, 1992; Nevin & Shahan, 2011) provides a quantitative framework for understanding why behavior persists in the face of disruption. The central finding is that behavior maintained by higher rates of reinforcement is more resistant to change, analogous to how a more massive object in motion is harder to stop. A key strength of the theory is that this same idea is expressed across different disruptors, each with its own equation.

Part 1: Prefeeding

You will work with data from a simulated multiple-schedule experiment in which subjects responded under both rich (60 reinforcers per hour) and lean (18 per hour) reinforcement components. After stable baselines were established, responding was disrupted via prefeeding at graded amounts (0, 25, 50, 75, and 100 g). Your task is to calculate proportion of baseline responding under each disruption level, fit the chapter's momentum equation with scipy.optimize.curve_fit, and test the theory's central claim: that one disruptor constant, combined with the known reinforcement rate in each component, accounts for the difference in resistance to change between the two components.

The equation you will fit is the chapter's:

log(Bx/B0) = -x * c / r

where Bx is the disrupted response rate, B0 is the baseline rate, x is the disruption level, c scales the effect of the disruptor, and r is the reinforcement rate in the component. Only c is unknown; x and r are set by the experimenter and are in the data file. You fit c separately for each component and then jointly with a single c, and compare.

Part 2: Other Disruptors

Prefeeding is only one way to disrupt responding. In this part you will fit the two further forms from the chapter (Nevin & Shahan, 2011; Nevin, McLean, & Grace, 2001) to data where the disruptor is extinction and where it is alternative reinforcement:

• Resistance to extinction: Bt/B0 = 10^(-t * (c + d*r) / sqrt(r)), with c and d the fitted constants.
• Alternative reinforcement: Bx/B0 = 10^(-p * Ra / sqrt(r + Ra)), with p the fitted constant.

Each data file holds several conditions with different, known baseline reinforcer rates r. From a single extinction curve only the slope (c + d*r)/sqrt(r) is identifiable, so you fit c and d jointly across conditions, passing both the disruptor and r to curve_fit as a two-row xdata array. The same applies to p.

Hints for Part 2:

• Take B0 as the session-0 (or Ra = 0) response rate in each condition and work with log10(Bt/B0).
• Constrain c, d, and p to be non-negative; start them at small positive values.
• Report standard errors from the covariance matrix curve_fit returns, and use sklearn.metrics for goodness of fit.

Through these exercises you will gain hands-on experience with the quantitative tools behavioral momentum theory offers for predicting and interpreting response persistence across qualitatively different disruptors.

Model Comparisons Lab

This week, we focus broadly on the topic of model comparisons. This week is less hands-on and more about ensuring that you understand the tools at your disposal. Models generally come in two flavors: regression (continuous or ordinal output) and classification (categorical output). The readings this week highlighted the many different ways we can think about comparing how close our model outputs are to the observed data (i.e., the various fit metrics).

When comparing models, we are choosing one (or more) fit metrics and seeing which model(s) lead to better fit metrics. The loss metric(s) we choose are often determined by our data, the audience, historical precedence, and potentially the downstream deployment environment. Unfortunately, there's no easy decision tree for making this decision. The "right" choice is often determined by thinking critically about your situation, what the model is doing, and how you can interpret the loss metric relative to the "meaning" of the data you are modeling. The nice thing is that choosing a loss metric is the hard part. Once a metric is chosen, comparing models is easy.

The only wrinkle is that models with a greater number of parameters have greater flexibility and, therefore, tend to have lower loss metric values. But the sacrifice is complexity and potentially transparency and explainability. This is where the information criteria metrics are useful (e.g., BIC, AIC, AICc). The idea with these is to penalize models for having more parameters to try to "even the scales." You should always include some measure of information criteria if you are comparing models with different numbers of parameters.

Part 1: Discounting Model Comparison (Regression)

The first dataset contains discounting data of the kind you fit in Week 2. This dataset will allow you to compare the following models of discounting: exponential, hyperbolic, Myerson and Green's hyperboloid, and Rachlin's hyperboloid. Fit each with scipy.optimize.curve_fit, as in Weeks 2, 3, and 5. For every row and for each of the four models, you should obtain the following fit metrics: r-squared, MAE, RMSE, AIC, BIC, and AICc. The first three come from sklearn.metrics; the information criteria are computed from the residual sum of squares with the chapter's formulas, because curve_fit does not report a likelihood.

Part 2: Challenging Behavior Prediction (Classification)

The second dataset is a held-out test set from a hypothetical classification task: predicting whether or not challenging behavior will occur during a session (yes=1, no=0). Three classifiers have already been fit to separate training data (a logistic regression, a decision tree, and a random forest), and the file gives you each model's predicted probability for every test session alongside the observed outcome. You will build these classifiers yourself in Week 12; this week the models are given and the question is how to compare them. Threshold the probabilities and obtain the following fit metrics from sklearn.metrics: accuracy, precision, recall, F1 score, Matthews Correlation Coefficient (MCC), and ROC-AUC. Then look at the confusion matrices and at how precision and recall trade off as the threshold moves. About 74% of sessions contain challenging behavior, so accuracy alone is not informative.

Useful Packages

scipy.optimize.curve_fit and sklearn.metrics (r2_score, mean_absolute_error, mean_squared_error, accuracy_score, precision_score, recall_score, f1_score, matthews_corrcoef, roc_auc_score, RocCurveDisplay, confusion_matrix).

How to Construct a Model Lab

This week, we focus on how to construct a model for some phenomenon you are interested in describing quantitatively. The chapters focused on how to build lifecycle diagrams and flow diagrams to start mapping out what variables are in your model and how they might relate. As with past chapters, we can then use the model we created to describe data we have collected, and use various techniques to improve those model fits (e.g., changing mathematically how the different variables relate).

From there, the chapters also covered how we might create a recursive model, a difference model, and a differential model, which allows us to run simulations to get a better feel for our models. This can be incredibly useful to determine when our model might make illogical predictions (e.g., response rates less than zero; exponentially growing response rates). It can also allow us to visualize how different parameters of our model influence our behavioral predictions.

During the in-class demo, we walked through two examples of how to do this for lever pressing under a random ratio (RR) schedule and salivary responding to a conditioned stimulus (CS). We also walked through how we can convert that equation into code to either fit data that we have collected or to run simulations. In each case, the models were one-variable models.

The notebook fits with scipy.optimize.curve_fit and scores fits with sklearn.metrics, the same tools as Weeks 3 and 5. One point the RR example makes explicit: when several parameters enter the model only through a single combination (there, $m = \alpha p - \beta$), data on the behavior alone cannot separate them, and the honest fit estimates the combination. Check your own models for the same issue before fitting. Next week you will refit the RR model by maximum likelihood, so keep your fitted values.

Assignment

Identify two different environment-behavior relations that you are interested in modeling. For each, create:

1. A life-cycle diagram
2. The basic equation relating environmental variables to behavior
3. The set of recursive, difference, and differential equations
4. The code to fit the model to a fake dataset
5. Simulations varying the parameters for each of the equations in step 3

Probability Theory and Probabilistic Models Lab

This lab works through two probabilistic tools from the chapter, each applied to behavioral data. Every estimate is obtained with a package function (statsmodels for the generalized linear model, scipy.stats for distributions, quantiles, and likelihoods); the chapter's hand derivations are there to show where those functions' answers come from, not to be repeated.

Part 1 takes the Poisson process from a rate estimate to a tested model. You will estimate a response rate by maximum likelihood with a Poisson GLM, read off its standard error and confidence interval, and then run the three checks the chapter describes: whether the variance of per-interval counts equals their mean, whether observed inter-response times follow the predicted exponential, and whether a quantile-quantile plot falls on the diagonal. Two 50-minute records are supplied that return the same rate estimate with overlapping confidence intervals. One was generated by a Poisson process and one by bouts of fast responding separated by pauses. Nothing in the point estimate distinguishes them; every diagnostic does. Part 1 closes with a bridge from Week 7: the random-ratio recursion you fit by least squares last week is a log-linear model in session number, so a Poisson GLM fits it by maximum likelihood in one call. You compare the two fits and state when the choice of likelihood matters.

Part 2 applies Bayesian updating to functional assessment. Starting from a uniform prior over four candidate functions (attention, escape, tangible, automatic), you will specify a Poisson likelihood encoding what each function predicts for each condition, then update after every session and track how the posterior moves. This is the chapter's four-hypothesis discrete example applied to real session records. Two datasets are supplied. In the clear case the posterior resolves within two sessions, and the sequence is worth following: a single elevated attention session leaves attention and automatic tied at 0.50 each, because both predicted elevation there, and the tie is broken by a condition in which responding was absent. In the ambiguous case the evidence per session is small, the posterior climbs gradually and non-monotonically, and a sensitivity analysis shows the final confidence depends on assumptions that the clear case was immune to.

Across both parts, a probabilistic model specifies a distribution rather than a single value, so the data can disagree with it in more than one way, and it states uncertainty as a quantity.

Packages

pandas, numpy, matplotlib, scipy.stats (poisson, expon, probplot), scipy.optimize.curve_fit, and statsmodels (smf.glm with sm.families.Poisson()).

Multilevel Modeling and Time-Series Forecasting Lab

Behavioral data are almost always nested: responses are nested within sessions, sessions are nested within participants, and participants may be nested within groups or settings. Ignoring this nesting structure when fitting statistical models can lead to biased parameter estimates and artificially narrow confidence intervals.

In this lab you will work with a dataset containing session-level response rates from eight participants, each observed across twenty sessions under varying reinforcement rates. You will begin by fitting a single-level ordinary least squares regression that ignores the nesting of sessions within participants. You will then fit a multilevel (mixed-effects) model that accounts for participant-level variation in both intercepts and slopes, compare the two approaches using information criteria, and interpret the fixed and random effects in behavioral terms.

You will then treat one participant's session-by-session rates as a time series, as in the chapter's Example B: detrend it, fit an AR(1) model with statsmodels.tsa.ar_model.AutoReg, check the residual autocorrelation, and forecast the next sessions. As an optional extension, you will decompose the same series into trend and residual components with seasonal_decompose.

Assignment

1. Load and explore the nested behavioral dataset.
2. Fit an OLS regression pooling all data and examine the residuals.
3. Fit a random-intercept mixed-effects model using statsmodels.formula.api.mixedlm.
4. Extend to a random-intercept-and-slope model.
5. Compare single-level and multilevel models using AIC and BIC.
6. Interpret the fixed effects and random effects in the context of the matching law and reinforcement rate-response rate relations.
7. Visualize participant-level regression lines.
8. Fit an AR(1) model to one participant's detrended series, interpret the autoregressive coefficient, and forecast three sessions ahead.
9. (Optional) Perform a simple time-series decomposition on one participant's data.

Dynamical Systems Models Lab

Many behavioral processes unfold over time in ways that are well described by differential equations. The logistic ordinary differential equation (ODE), dx/dt = r * x * (1 - x/K), provides a simple but powerful model of acquisition: responding starts slowly, accelerates as the behavior contacts reinforcement, and then decelerates as it approaches a carrying-capacity asymptote.

In this lab you will implement the logistic ODE numerically using Euler's method (once, so you know what a solver does), then solve the same equation with scipy.integrate.solve_ivp, the adaptive solver you would use in practice, and compare both against the closed-form solution. You will explore how the parameters r (growth rate) and K (carrying capacity) shape the acquisition curve, find the equilibrium points analytically, perform a linear stability analysis, and construct a phase portrait. You will then fit the analytical solution of the logistic equation to empirical acquisition data with scipy.optimize.curve_fit, scoring the fit with sklearn.metrics.

Assignment

1. Implement Euler's method for the logistic ODE and plot the resulting acquisition curves.
2. Vary the initial condition x0 and the parameters r and K to see how they affect the shape and speed of acquisition.
3. Find the equilibrium points of the logistic ODE analytically.
4. Perform a linear stability analysis by evaluating f'(x*) at each equilibrium.
5. Create a phase portrait (dx/dt vs. x) with directional arrows.
6. Solve the ODE with scipy.integrate.solve_ivp and compare its error and cost against Euler's method.
7. Load the empirical acquisition dataset and fit the logistic analytical solution using scipy.optimize.curve_fit.
8. Plot the fitted curve alongside the data and report the estimated parameters.
9. Discuss the biological and behavioral interpretation of r and K.

Computational Models Lab

This week moves from descriptive models of behavior to process models that specify how behavior is generated. Rather than fitting equations to existing data, you will build an agent that learns through trial-and-error interaction with a simulated concurrent VI 30-s VI 60-s schedule.

The lab repeats the chapter's worked example. You will implement the schedule (one-second trials, a reinforcer set up on each lever with probability 1/VI and held until collected) and a Q-learning agent with softmax action selection, using the chapter's parameters ($\alpha = 0.1$, $\gamma = 0.95$, $\tau = 0.1$). Each session is 3,600 trials with a new agent, and steady state is measured over the second half. You will run 100 sessions, compare the agent's allocation with the matching prediction, and use the chapter's fixed-point result to explain why it undermatches.

You will then sweep the learning rate, the discount factor, and the softmax temperature to see which parameters move steady-state allocation, and fit the generalized matching equation from Week 2 to the agent's allocation across five schedule pairs with scipy.stats.linregress.

The last two tasks apply the chapter's Monte Carlo material to your own simulation. A single run returns a single estimate of sensitivity, and that estimate has sampling variability that one run hides. You will replicate the whole five-pair procedure 100 times to obtain the Monte Carlo distribution of sensitivity. Then you will treat 10 sessions as if they were the only data, use scipy.stats.bootstrap to put a confidence interval on the mean steady-state choice proportion, and compare it with the Week 8 interval, the estimate $\pm 1.96 \times SE$. The bootstrap is the tool for the case where you cannot rerun the process, which is every empirical data set.

The environment and agent are written from scratch because the update rule is what this week is about. They interact only through reset() and step(), the common interface pattern for reinforcement learning environments.

Assignment

Complete all tasks in the Jupyter notebook. Each task builds on the previous one. Pay particular attention to the final discussion, where you reflect on the relationship between descriptive models (matching law) and process models (Q-learning).

Machine Learning and Artificial Intelligence Lab

This week introduces supervised machine learning as a tool for prediction and classification in behavior science. The specific application is classifying the function of problem behavior from functional analysis (FA) summary data -- a task that clinicians perform routinely based on visual inspection of FA graphs.

You will work with a simulated dataset of 120 participants, each characterized by response rates across four FA conditions (attention, escape, tangible, and play/control). Each participant's behavior is maintained by attention, escape, or tangible reinforcement, and the groups overlap enough that no classifier is perfect. Your goal is to build a classifier that can predict the behavioral function from these features. You will start with a single decision tree, which is highly interpretable and mirrors the kind of rule-based reasoning clinicians use (e.g., "if the rate in the attention condition is elevated relative to play, the function is likely social-positive reinforcement").

From there, you will explore overfitting by comparing training accuracy with cross-validated accuracy across different tree depths (choosing the depth by cross-validation on the training set, as in Week 6, rather than by peeking at the test set), then build a random forest to see whether an ensemble of trees improves predictive performance. You will then fit a logistic regression as the linear baseline and compare all three classifiers with the metrics you used on pre-fit models in Week 6 (accuracy, macro F1, MCC, and per-class precision and recall). The lab concludes with a discussion of the prediction-explanation tradeoff: more complex models may predict better, but simpler models are easier to interpret and communicate to practitioners.

Assignment

Complete all tasks in the Jupyter notebook. You will need scikit-learn, pandas, and matplotlib. Focus not only on model accuracy but on understanding why the models make the predictions they do.