6.4 经典漂移扩散模型

Edited by Yuhong

前一小节中，我们通过假设决策过程中存在漂移-扩散的过程，模拟了随机点运动实验中决策边界、正确率和反应时的关系。我们可以发现，速度-准确性的权衡受到决策边界的影响：越高的决策边界会导致越长的反应时和正确率，而越低的决策边界会导致更短的反应时和更差的正确率。这意味着个体决策边界的不同会影响数据的解释，然而在数据分析中往往忽略了这一问题。解决这个问题的一种方法是使用计算模型来提取潜在的认知过程的参数。在本章中，我们将介绍提取决策边界参数的其中一种模型——漂移-扩散模型，并以该思想在python中拟合随机点运动的实验。但在这之前，让我们重新回顾一下心理学计算模型参数与认知过程的关系。

由于认知过程往往不能被直接观测，认知心理学家们提出将人脑类比于计算机，提出了信息加工理论来解释人类底层的认知过程。信息加工理论认为，认知过程像计算机一样，按时间顺序固定阶段接受、存储、处理并输出信息。例如，在玩你画我猜的游戏时，猜图像的人首先需要对图像进行感知，然后通过注意进入工作记忆系统，进而进行进一步的加工、决策，最后转化为动作（语言）输出。受到这种思想的启发，心理学家们经常通过外显的行为或神经测量数据来构建认知计算模型，以解释我们的认知过程，例如“决策边界”和“信息累积速度”。其中，认知计算模型参数与这些认知机制相关联。我们可以把真实数据的分布（例如，倒过来的“U”型分布）想象为在玩“你画我猜”游戏时，出题者所欲让别人猜出来的图形。认知计算模型可以理解为出题者实际画出来的图形，而模型参数则决定我们画出来的图形长什么样。调整模型参数时，我们画出来的的图形和真实图形会更相似或者更不相似。当我们修改参数，使模型的图像尽可能接近实际数据，就是“将模型与数据拟合”的过程。这样做的好处是可以将实际数据与认知理论中的抽象概念建立起更加直观的联系，从而帮助我们理解和预测人类的认知过程。通过调整模型参数来拟合数据，我们不仅可以以通过模型比较，找出最佳模型验证已有的理论假设，并且还能量化隐藏在数据背后的潜在认知机制。

漂移-扩散模型(drift-diffusion model, DDM)是目前在知觉决策领域最广为使用的模型(Ratcliff, 1978)。自从被Ratcliff于上个世纪提出，经过60多年的发展，漂移-扩散模型已在多个领域被证明有良好的有效性，包括记忆、知觉决策、社会认知和精神疾病等，并且具有良好的神经证据(Gold & Shadlen, 2007)。漂移-扩散模型是证据积累模型（又名序列采样模型）的一种，它将决策假设为随着时间推移对来自感觉或者记忆的刺激连续不断地积累证据的过程，当支持一个反应的累积证据超过阈值时做出反应。在经典的漂移扩散模型中，反应时间被拆解为非决策时间($T_{er}$, encode and response time)和决策时间($T_d$, decision time)。其中非决策时间可以拆分为刺激编码($T_e$, encode time)和动作反应($T_r$, response time)的时间，决策时间则代表证据不断积累直到到达其中一个反应的边界的时间。在决策时间（在刺激被编码之后，但在做出决定的反应之前）中，证据从起始点($z$, initial bias)开始以带有噪声的方式积累。这种噪声来源于在决策过程中每一瞬间对特定决策选项的新的证据积累，表现形式为每一时刻累积证据都会随机向上或下边界游走。这种噪声被认为模拟了决策过程中大脑神经元放电的随机性和注意力的波动(Myers et al., 2022)。虽然证据的积累过程带有噪声，但最终会到达上或下边界($a$, boundary separation)以做出决策。其中上边界代表正确的响应，范围是(0, +∞)；下边界代表错误的响应，一般为0。在这一段证据积累的过程中，证据积累的平均速度是漂移率($v$, drift rate)。

（图1）

由此，我们对经典漂移-扩散模型的4个参数进行总结：

参数名称	参数缩写	值的范围	代表的认知过程
漂移率（Drift Rate）	v或 d(Myers et al., 2022)）	(0, +∞)	证据积累的速度，可以用来反映任务的难度之间的差异
决策边界 （Boundary Separation）	a	(0, +∞)	决策所需要的证据量，反映了决策者的速度-准确性权衡。决策阈限越大意味着越慢但越准确的响应，决策阈限越小意味着越快但错误率越高的响应。
决策起点 （Initial Bias）	z	(0, a)	响应偏差，代表在信息开始积累之前对其中一个反应的偏好。当z = a /2时代表反应是无偏的，当z < a /2时代表反应偏向下边界，当z > a /2时代表反应偏向上边界。
非决策时间 (No Decision Time)	t或 Ter(Weigard et al., 2021)、$T_er$(Ratcliff et al., 2016)、$t_0$ (Voss et al., 2004)）	(0, +∞) （不可以超过总反应时）	决策过程之外的时间，代表对刺激进行感知编码和做出动作反应的过程。

接下来，让我们使用经典的漂移扩散模型思想，再次模拟随机点运动的实验。在进行模拟之前，我们需要对模型的参数进行申明，并介绍拟合模型的方法。

首先，在标准的漂移-扩散模型中，假设每一刻所积累的证据满足正态分布。我们假设在随机点运动实验中，每一帧所积累的感官证据的满足的以上正态分布的宽度（即标准差）σ = 0.1。由于朝同一个方向运动的点的一致性百分比(coherence)越高，刺激越明显，每一帧感官证据积累的强度越强，所以整体的证据积累会越快，所以我们假定每一帧积累的感官证据的正态分布均值$\mu$和coherence为一个线性关系：$\mu=k*coherence$。其中，$k$是漂移系数，代表个体信息累计的能力，$\mu$代表证据累计速度。与前几节的设定不同，我们按经典的漂移扩散模型设定两个决策边界，其中上边界为正确的响应，下边界(y = 0)为错误的响应。假设证据从决策边界的一半开始，从初始点到上边界和下边界的距离是一样的。此外，我们假设在随机点运动中有N个试次，我们已知的数据有每个试次的：

1. Coherence

2. 正确还是错误(correct/wrong, 1/0)

3. 反应时 reaction time (milliseconds)

但我们不知道被试的四个漂移-扩散模型参数（见表x），所以我们根据数据来反推四个参数。我们令$θ = (k,B,z,t)$，其中$θ$是参数向量，代表我们所不知道的参数；k是漂移系数，在此根据我们之前作出的假设替代了漂移率v；a是决策边界；t是非决策时间。

在确定了未知的参数和数据后，接下来我们要将模型与数据拟合，以估计具体的参数值。目前有多种成熟的方法可以进行参数估计，如：$\Chi^2$卡方检验(Ratcliff & Tuerlinckx, 2002)，最大似然估计（见第三章第二小节），贝叶斯方法(Kruschke & Liddell, 2018)和贝叶斯分层方法(Vandekerckhove et al., 2011)。每种方法都有其各自的优点和缺点，我们推荐使用贝叶斯分层方法，但因为使用python代码进行复现较为困难，在演示中我们使用最大似然估计方法进行拟合：

从概率推断的角度，我们用最大似然估计的方法，假定有N个试次，已知第i个试次的coherence 水平$coh_i$，每个试次的选择是正确还是错误$c_i$，以及反应时$rt_i$。对于所有试次来说，负对数似然函数（见本书3.2.2(3.2)）是：$\Chi^2$

$$\text{NLL} = -LL(\theta) = -\sum_{i=1}^{N} \log p(rt_i \mid \theta, \ coh_i, \ c_i)$$

其中：

$\text{NLL}= -LL(\theta)$是负对数似然函数。我们希望通过调整参数$θ$来最小化这个损失函数，从而找到最符合数据的参数组合。

$p(rt_i \mid \theta, \ coh_i, \ c_i)$指的是在给定参数$θ$、一致性$coh_i$、和决策$c_i$的情况下，模型预测在第i次试验中观测到反应时间$rt_i$的概率。

$\sum_{i=1}^{N} \log$表示我们要对所有的N个试次的观察值计算似然函数。逐个试次计算反应时间和决策的概率，然后取对数累加。

首先，让我们用漂移-扩散模型思想模拟一个试次中的证据积累过程。

我们先对部分参数做出定义：

```python
from scipy.stats import norm #用于从正态分布中采样随机数。这里正态分布代表每一帧输入的证据波动。
import numpy as np #用来做数值计算，比如生成数组、计算随机数等。
import matplotlib.pyplot as plt #用来绘制图表，帮助我们可视化累积的证据变化过程。

B = 2 # 我们假设底部的边界为0，顶部的边界为2
k = 1 # 漂移率和朝同一个方向运动的点的一致性百分比之间的比例系数
coh = 0.064 # choose one coefficient #朝同一个方向运动的点的一致性百分比，越低任务越难
DriftRate = k * coh # 漂移率
sigma = 0.1 #正态分布的宽度（即标准差）

然后，让我们对数据进行初始化：

evidence = B/2 #设定初始证据为一半的边界值，假设一开始对正确和错误的决策倾向是无偏的
accumEvidence = [] #创建列表以储存每一帧积累的证据
accumEvidence.append(evidence)

接着，让我们循环模拟决策过程中的证据积累：

```python
i = 0 #记录当前是第几帧
while np.abs(evidence) < B: #当累积的证据绝对值还没达到决策边界a时，继续累积证据，直到累积到边界停止
    evidence = evidence + norm.rvs(loc=driftRate, scale=sigma, size=1)[0] #从正态分布中生成一个随机数，用来模拟漂移-扩散模型中的随机游走
    accumEvidence.append(evidence)
    i = i+1 

最后，让我们来对一个试次中的证据积累过程进行可视化：

```python
plt.figure(figsize=(4,4))
plt.plot(np.arange(i+1), accumEvidence, '-o')
plt.axhline(B, color='k', label='Decision boundary')
plt.xlabel('# of frame')
plt.ylabel('Accumulated evidence')
plt.ylim([0, B+1])

在用漂移-扩散模型模拟完一个试次中的证据累计过程后，让我们接着使用最大似然估计值，根据真实的数据来估计参数值。首先，让我们先构建经典漂移-扩散模型中的概率密度函数：

```python
from math import pi, sqrt, log, ceil, floor, exp, sin
from numpy import log


def ftt_01w(tt, w, err=1e-4):
    """根据 Navarro 和 Fuss (2009) 计算 f(t|0,1,w)"""
    # tt 标准化的反应时间
    # w 起始点，以决策边界的比例代替
    # err 允许的计算误差

    # 计算大t所需的项数k1，表示需要多少项才能满足误差要求
    if pi * tt * err < 1:  # 检查错误阈值是否足够低，以使项数近似足够准确。
        kl = sqrt(-2 * log(pi * tt * err) / (pi ** 2 * tt))  # 使用对数近似来估计误差足够小时所需的项数
        kl = max(kl, 1.0 / (pi * sqrt(tt)))  # 确保满足准确度所需的最小项数
    else:  # 如果错误阈值设定的太高
        kl = 1.0 / (pi * sqrt(tt))  # 设定边界条件

    # 计算小 t 所需的项数，表示需要多少项才能满足误差要求
    if 2 * sqrt(2 * pi * tt) * err < 1:  # 检查错误阈值是否足够低，以使项数近似足够准确。
        ks = 2 + sqrt(-2 * tt * log(2 * sqrt(2 * pi * tt) * err))  # 使用对数近似来估计误差足够小时所需的项数
        ks = max(ks, sqrt(tt) + 1)  # 确保满足准确度所需的最小项数
    else:  # 如果错误阈值设定的太高
        ks = 2  # 默认使用2个项

    # 计算 f(tt|0,1,w)
    p = 0  # 初始化密度
    if ks < kl:  # 如果小t更好
        K = ceil(ks)  # 将ks向上取整，表示我们在循环中需要处理的最大迭代次数
        lower = -floor((K - 1) / 2.0) # 计算循环的上界
        upper = ceil((K - 1) / 2.0)  # 计算循环的下界
        for k in range(lower, upper + 1):  # 在从lower到upper界限之间的每一个整k上循环，计算每个k对p的贡献
            p += (w + 2 * k) * exp(-(pow((w + 2 * k), 2)) / 2 / tt)  # 对p进行累加
        p /= sqrt(2 * pi * pow(tt, 3))  # 添加常数项

    else:  # 如果大t更好
        K = ceil(kl)  # 将kl向上取整，表示我们在循环中需要处理的最大迭代次数
        for k in range(1, K + 1):  # 大时间情况下，从1到k的循环
            p += (
                    k * exp(-(pow(k, 2)) * (pi ** 2) * tt / 2) * sin(k * pi * w)
            )  # 累加每个k对p的影响
        p *= pi  # 添加常数项

    return p


def pdf(x, k, B, a, err=1e-4):
    """根据 Navarro 和 Fuss (2009) 计算 f(x|k,B,a)"""
    # 反应时间必须是正数
    # x 反应时间
    # k 漂移系数
    # B 决策边界
    # a 起始点，以决策边界的比例代替
    # err 允许的计算误差

    if x <= 0:  # 如果反应时间x小于或等于0，函数直接返回0，因为反应时间不可能为负数
        return 0

    tt = x / B ** 2  # 将反应时间标准化

    p = ftt_01w(tt, w=a, err=err)  # 调用ftt_01w函数，使用标准化时间tt和偏差a来计算概率密度p

    # 转化为f(t|v,a,w)
    xx = exp(-k * B * a - (pow(k, 2)) * x / 2.0) / (pow(B, 2))  # 计算转化系数xx，用来将标准化的结果转换为有漂移情况下的
    r = p * xx  # 将概率密度p乘以转换系数xx，计算得到最终的概率密度值r
    # print(x, v, a, w,p,r,xx,tt)
    return r

def ddmpdf(k, a, B, ndt, coh, correct, rt):
    """
    该函数接受七个参数, 其中四个是ddm的参数
    k: 漂移系数
    a: 初始偏差, 这里表示为和B的比例, 范围是(0,1)
    B: 决策边界
    ndt: 非决策时间

    另外三个是一个trial的数据
    coh: coherence
    correct: 是否为正确的决策 (1 正确, 0 错误)
    rt: 反应时间，以秒为单位
    """
    if correct==0: # 错误的trial
        return pdf(x=rt-ndt, k=k*coh, B=B, a=a) #调用pdf函数，使用经过非决策时间ndt修正的反应时间，漂移系数k根据coh进行调整。初始偏差a保持不变
    elif correct==1: # 正确的trial
        return pdf(x=rt-ndt, k=-k*coh, B=B, a=1-a) # 也调用pdf函数，但漂系数k取负值，同时初始偏差a改为1-a（反向累积）
```

在定义好经典漂移-扩散模型的概率密度函数后，让我们先导入示例数据：

```python
import requests
import numpy as np
from io import StringIO

# 从 URL 获取数据
response = requests.get('https://raw.githubusercontent.com/sunyuhongwr/-/main/DDMexampledata.txt')
# 使用 StringIO 处理响应内容
data = np.loadtxt(StringIO(response.text), delimiter="\t")
# 打印数据的整体情况
print(data.shape)
```

(300, 3)

实例数据是一个300（试次）x 3列 的数组。其中：

第一列：coherence

第二列：反应时间

第三列：正确/错误（1/0）

接着，让我们定义负对数似然函数，用最大似然估计法估计参数：

```python
#  定义负对数似然函数
def negloglikeli(params):
    '''
    <params包含的参数>: 漂移系数, 决策边界, 初始偏差, 非决策时间
    '''
    k = params[0] # 漂移系数    
    B = params[1] # 决策边界
    a = params[2] # 初始偏差 (0, 1)
    ndt = params[3] # 非决策时间

    nTrial = data.shape[0] # 获取数据中的试次数量
    pp = np.empty(nTrial)  # 初始化一个空的数组pp,用于存储每个试次的概率密度值
    for i in range(nTrial): # 循环每一个试次
        pp[i] = ddmpdf(k, a, B, ndt, coh=data[i, 0], rt=data[i, 1], correct=data[i, 2])  #通过ddmpdf函数计算每个试次的概率密度值pp[i]

    pp=0.999*pp + np.finfo(np.float32).eps # 避免概率为0
    return -np.log(pp).sum()  # 取对数并累加，然后返回负值
```

最后，让我们优化这个函数来获得拟合的参数

 ```python
from scipy.optimize import minimize
from ddm import ddmpdf # 漂移-扩散模型的概率密度函数，我们先前已作出定义
res = minimize(fun=negloglikeli, x0=(1, 2, 0.5, 0.2), bounds=((0, 10), (0, 5), (0, 1), (0, 1))) # 使用 scipy.optimize.minimize 函数进行参数优化（使负对数似然函数最小化的漂移系数、决策边界、初始偏差和非决策时间）；x0指初始参数值，其中四个数值分别代表漂移系数、决策边界、初始偏差和非决策时间；bounds定义了边界来限制参数的范围，其中四个数值分别代表漂移系数、决策边界、初始偏差和非决策时间

print('\n拟合的漂移系数是 ', res.x[0])
print('拟合的决策边界是 ', res.x[1])
print('拟合的起始偏差是 ', res.x[2])
print('拟合的非决策时间是 ', res.x[3])
```

拟合的漂移系数是 5.230241590116383

拟合的决策边界是 3.7366640686792887

拟合的起始偏差是 0.5365699686795544

拟合的非决策时间是 0.03953784720249189

虽然在本次演示中的代码较为复杂，但是在实际使用漂移-扩散模型中已有较为成熟的工具包来计算。具体有以下一些：

Shinn, M., Lam, N. H., & Murray, J. D. (2020). A flexible framework for simulating and fitting generalized drift-diffusion models. eLife, 9, e56938. https://doi.org/10.7554/eLife.56938

参考文献

Gold, J. I., & Shadlen, M. N. (2007). The Neural Basis of Decision Making. Annual Review of Neuroscience, 30(1), 535–574. https://doi.org/10.1146/annurev.neuro.29.051605.113038

Kruschke, J. K., & Liddell, T. M. (2018). Bayesian data analysis for newcomers. Psychonomic Bulletin & Review, 25(1), 155–177. https://doi.org/10.3758/s13423-017-1272-1

Myers, C. E., Interian, A., & Moustafa, A. A. (2022). A practical introduction to using the drift diffusion model of decision-making in cognitive psychology, neuroscience, and health sciences. Frontiers in Psychology, 13, 1039172. https://doi.org/10.3389/fpsyg.2022.1039172

Ratcliff, R. (1978). A theory of memory retrieval. Psychological Review, 85(2), 59–108. https://doi.org/10.1037/0033-295X.85.2.59

Ratcliff, R., Smith, P. L., Brown, S. D., & McKoon, G. (2016). Diffusion Decision Model: Current Issues and History. Trends in Cognitive Sciences, 20(4), 260–281. https://doi.org/10.1016/j.tics.2016.01.007

Ratcliff, R., & Tuerlinckx, F. (2002). Estimating parameters of the diffusion model: Approaches to dealing with contaminant reaction times and parameter variability. Psychonomic Bulletin & Review, 9(3), 438–481. https://doi.org/10.3758/BF03196302

Vandekerckhove, J., Tuerlinckx, F., & Lee, M. D. (2011). Hierarchical diffusion models for two-choice response times. Psychological Methods, 16(1), 44–62. https://doi.org/10.1037/a0021765

Voss, A., Rothermund, K., & Voss, J. (2004). Interpreting the parameters of the diffusion model: An empirical validation. Memory & Cognition, 32(7), 1206–1220. https://doi.org/10.3758/BF03196893

Weigard, A., Clark, D. A., & Sripada, C. (2021). Cognitive efficiency beats top-down control as a reliable individual difference dimension relevant to self-control. Cognition, 215, 104818. https://doi.org/10.1016/j.cognition.2021.104818