{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Behavioral Momentum and Response Persistence\n",
    "\n",
    "In this lab, you will analyze resistance to change across reinforcement contexts, using three different *disruptors*. Behavioral momentum theory (Nevin, 1992; Nevin & Shahan, 2011) predicts that responding maintained by richer reinforcement is more resistant to disruption.\n",
    "\n",
    "The lab has two parts:\n",
    "\n",
    "- **Part 1 - Prefeeding.** Fit the simplified momentum equation to multiple-schedule data where prefeeding is the disruptor, and compare resistance to change across rich and lean components.\n",
    "- **Part 2 - Other disruptors.** Fit Nevin's resistance-to-extinction and alternative-reinforcement equations to data where extinction and alternative reinforcement are the disruptors.\n",
    "\n",
    "## Background (Part 1)\n",
    "\n",
    "Four subjects responded on a two-component multiple schedule with **rich** and **lean** reinforcement components. After baseline response rates were established, responding was disrupted by prefeeding at graded levels (0, 25, 50, 75, and 100g). The data record baseline and disrupted rates for each subject, component, and disruption level.\n",
    "\n",
    "The behavioral momentum equation describes the relation between disruption and proportional change in responding:\n",
    "\n",
    "$$\\log\\left(\\frac{B_x}{B_0}\\right) = \\frac{-x \\cdot c}{r \\cdot S}$$\n",
    "\n",
    "where $B_x$ is the response rate under disruption level $x$, $B_0$ is the baseline rate, $c$ is sensitivity to disruption, $r$ is the reinforcement rate, and $S$ captures the stimulus-reinforcer relation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 1: Import Libraries\n",
    "\n",
    "Import the libraries you will need for this lab: `pandas`, `numpy`, `matplotlib.pyplot`, and `scipy.optimize`."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Import your libraries here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 2: Load and Inspect the Data\n",
    "\n",
    "Load `momentum_data.csv` into a DataFrame. Examine the structure of the data. How many subjects are there? What are the unique disruption levels? Print summary statistics for the rich and lean components separately."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Load momentum_data.csv and inspect it\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 3: Calculate Proportion of Baseline\n",
    "\n",
    "For each row in the dataset, calculate the proportion of baseline responding ($B_x / B_0$). Add this as a new column called `prop_baseline`. Then calculate $\\log_{10}(B_x / B_0)$ and store it in a column called `log_prop_baseline`.\n",
    "\n",
    "**Note:** At disruption level 0, the proportion should be 1.0 and the log proportion should be 0.0. Verify this is the case in your data."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Compute prop_baseline and log_prop_baseline\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 4: Visualize the Raw Disruption Functions\n",
    "\n",
    "Create a figure with two panels (subplots). In the left panel, plot proportion of baseline ($B_x / B_0$) as a function of disruption level, with separate lines for each subject. Use different colors or markers for rich vs. lean components. In the right panel, plot the same data on a log scale ($\\log_{10}(B_x / B_0)$ on the y-axis).\n",
    "\n",
    "Add appropriate axis labels, legends, and a title. What pattern do you notice regarding the rich vs. lean components?"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Two-panel plot: prop_baseline and log_prop_baseline vs disruption level\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 5: Compute Group Means\n",
    "\n",
    "Calculate the mean proportion of baseline and mean log proportion of baseline across subjects for each combination of component (rich/lean) and disruption level. Store these in a new DataFrame. This will be useful for fitting the momentum equation to averaged data."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Group means by component and disruption level\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 6: Define the Behavioral Momentum Equation\n",
    "\n",
    "Write a Python function that implements the behavioral momentum equation:\n",
    "\n",
    "$$\\log\\left(\\frac{B_x}{B_0}\\right) = \\frac{-x \\cdot c}{r \\cdot S}$$\n",
    "\n",
    "The function should take `x` (disruption level) as the independent variable and `c`, `r`, and `S` as parameters. It should return the predicted log proportion of baseline.\n",
    "\n",
    "**Hint:** Since $r$ and $S$ always appear as the product $r \\cdot S$, you may find it easier to define a combined parameter (e.g., `rS`) to avoid identifiability issues during fitting. Think about what this means for interpreting your results."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Define momentum_model(x, c, rS)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 7: Fit the Model to Each Component\n",
    "\n",
    "Using `scipy.optimize.curve_fit`, fit the behavioral momentum equation to the group-mean data for the **rich** component and the **lean** component separately.\n",
    "\n",
    "Report the best-fitting parameter values for each component. Which component shows a steeper disruption function (i.e., less resistance to change)?\n",
    "\n",
    "**Hints:**\n",
    "- Exclude the disruption level 0 data point from fitting (it is trivially 0 on the log scale).\n",
    "- Provide reasonable initial parameter guesses.\n",
    "- Consider whether parameter bounds are needed."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Fit momentum_model to rich and lean group means (exclude level 0)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 8: Evaluate Model Fit\n",
    "\n",
    "For each component, calculate the following goodness-of-fit measures:\n",
    "\n",
    "1. **R-squared ($R^2$):** proportion of variance accounted for.\n",
    "2. **Root Mean Squared Error (RMSE):** average magnitude of the residuals.\n",
    "\n",
    "How well does the behavioral momentum equation describe the disruption data? Are the fits comparable across components?"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Compute R-squared and RMSE per component\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 9: Plot Model Predictions Against Data\n",
    "\n",
    "Create a publication-quality figure showing:\n",
    "- The group-mean log proportion of baseline data points for each component (use distinct markers for rich and lean).\n",
    "- The fitted model predictions as smooth curves overlaid on the data.\n",
    "- A legend identifying the components.\n",
    "- Appropriate axis labels (\"Disruption Level (g prefeeding)\" and \"log(Bx/B0)\").\n",
    "\n",
    "Generate the smooth curves by evaluating the model function at many x-values between 0 and 100."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Plot data points and fitted curves for rich and lean\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 10: Interpretation and Discussion (Part 1)\n",
    "\n",
    "In a markdown cell below, address the following:\n",
    "\n",
    "1. How do the fitted parameter values differ between the rich and lean components? What do these differences mean in terms of behavioral momentum theory?\n",
    "2. Is resistance to change greater in the rich component, as predicted by the theory? How do you quantify this from your model fits?\n",
    "3. What are the limitations of using the simplified momentum equation with prefeeding as the disruptor?\n",
    "4. If you combined the parameter $r \\cdot S$ into a single value, what additional information would you need to separate them? Why might this matter for applied predictions?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2: Other Disruptors\n",
    "\n",
    "### Task 11: Resistance to Extinction\n",
    "\n",
    "Prefeeding is only one way to disrupt responding. Load `behavioral_momentum_extinction_data.csv`, where **extinction** is the disruptor across sessions. Implement Nevin's resistance-to-extinction equation\n",
    "\n",
    "$$\\frac{B_t}{B_0} = 10^{\\,-t\\,(c + d\\,r)/\\sqrt{r}}$$\n",
    "\n",
    "and use `curve_fit` to estimate `c`, `d`, `r`, and `B0` for each condition. Report the fitted parameters and $R^2$, and overlay the fitted curves on the data. (Bounds: $c, d > 0$; $r > 0.1$; $B_0$ a plausible baseline rate.)"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Load extinction data; define momentum_extinction(t, c, d, r, B0); fit each condition\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 12: Disruption by Alternative Reinforcement\n",
    "\n",
    "Load `behavioral_momentum_alternative_data.csv`, where responding is disrupted by **alternative reinforcement** delivered at rate $R_a$. Implement\n",
    "\n",
    "$$\\frac{B_x}{B_0} = 10^{\\,-p\\,R_a/\\sqrt{r + R_a}}$$\n",
    "\n",
    "and fit `p`, `r`, and `B0` per condition. Overlay the fitted curves, and quantify how much responding is suppressed from baseline to the highest alternative-reinforcement rate."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Load alternative-reinforcement data; define momentum_alternative(Ra, p, r, B0); fit each condition\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Wrap-up\n",
    "\n",
    "Across all three disruptors (prefeeding, extinction, alternative reinforcement), what is the common signature of behavioral momentum? Compare how each disruptor enters its equation. Which parameters index *sensitivity to disruption* versus the *baseline mass* of behavior, and how do the recovered values line up with the theory's prediction that richer reinforcement yields greater resistance to change?"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}