NeuromatchAcademy课程：贝叶斯决策模型拟合实战指南

前言

本教程是NeuromatchAcademy课程"贝叶斯决策"系列的第三部分，将深入探讨如何将贝叶斯模型拟合到实验数据中。我们将通过数学推导和代码实现，完整展示从模型构建到参数估计的整个过程。

模型概述

我们研究的场景是：参与者需要估计隐藏在幕布后的木偶发出的声音位置。参与者被告知声音可能来自两种原因：

1. 共同原因（95%概率）：声音来自幕布后位置0处的木偶
2. 独立原因（5%概率）：声音来自远处的扬声器

这个实验设置与第二教程类似，但概率参数有所不同。我们将构建一个包含高斯先验和高斯似然的贝叶斯混合模型。

核心概念

生成模型流程

1. 刺激呈现：向参与者呈现真实刺激x
2. 神经编码：大脑对刺激进行噪声编码，形成神经表征p(x̃|x)
3. 贝叶斯整合：将神经表征(似然)与先验信息p(x)结合，形成后验p(x|x̃)
4. 行为响应：基于后验估计产生行为响应x̂

关键步骤

我们将分两步实现模型拟合：

1. 构建模型组件：先验、似然、后验等
2. 执行模型反演/拟合：估计生成数据的模型参数

实现细节

1. 构建似然阵列

首先我们需要创建似然函数，但为了可视化和考虑所有可能的神经编码，我们将为每个可能的编码刺激x̃创建一个似然函数p(x̃|x)。

import numpy as np
import matplotlib.pyplot as plt

def my_gaussian(x_points, mu, sigma):
    """计算高斯分布"""
    p = np.exp(-(x_points-mu)**2/(2*sigma**2))
    return p / sum(p)

def compute_likelihood_array(x_points, stim_array, sigma=1.):
    """计算似然阵列"""
    likelihood_array = np.zeros((len(stim_array), len(x_points)))
    for i in range(len(stim_array)):
        likelihood_array[i, :] = my_gaussian(x_points, stim_array[i], sigma)
    return likelihood_array

# 定义参数
x = np.arange(-10, 10, 0.1)
hypothetical_stim = np.linspace(-8, 8, 1000)

# 计算并绘制似然阵列
likelihood_array = compute_likelihood_array(x, hypothetical_stim)

这段代码会生成一个2D阵列，其中每行是对应不同均值的高斯分布，展示了在不同真实刺激x下，神经编码x̃的概率分布。

2. 构建先验分布

接下来我们实现混合先验分布，包含两个成分：

• 共同原因：集中在0位置的高斯分布
• 独立原因：均匀分布

def compute_prior_array(x_points, p_common, prior_sigma=1.):
    """计算先验分布"""
    prior_common = my_gaussian(x_points, 0, prior_sigma)
    prior_independent = np.ones_like(x_points) / len(x_points)
    prior_array = p_common * prior_common + (1 - p_common) * prior_independent
    return prior_array

3. 计算后验分布

根据贝叶斯定理，后验分布是似然和先验的乘积：

def compute_posterior_array(likelihood_array, prior_array):
    """计算后验分布"""
    posterior_array = np.zeros_like(likelihood_array)
    for i in range(likelihood_array.shape[0]):
        posterior_array[i, :] = likelihood_array[i, :] * prior_array
        posterior_array[i, :] /= np.sum(posterior_array[i, :])
    return posterior_array

4. 行为响应模型

我们假设参与者的响应x̂是后验分布的均值：

def compute_resp(posterior_array, x_points):
    """计算响应分布"""
    resp_array = np.zeros(posterior_array.shape[0])
    for i in range(posterior_array.shape[0]):
        resp_array[i] = np.sum(x_points * posterior_array[i, :])
    return resp_array

模型拟合实战

数据生成

首先我们使用生成模型创建模拟数据：

def generate_data(n_samples, p_common, sigma_likelihood=1.):
    """生成模拟数据"""
    true_stim = np.random.uniform(-8, 8, n_samples)
    behaviour = np.zeros_like(true_stim)
    
    for i in range(n_samples):
        # 神经编码
        x_tilde = np.random.normal(true_stim[i], sigma_likelihood)
        
        # 贝叶斯估计
        prior = compute_prior_array(x, p_common)
        likelihood = my_gaussian(x, x_tilde, sigma_likelihood)
        posterior = likelihood * prior
        posterior /= np.sum(posterior)
        
        # 行为响应
        behaviour[i] = np.sum(x * posterior)
    
    return true_stim, behaviour

模型反演

现在我们进行模型反演，即从行为数据中估计p_common参数：

def neg_log_likelihood(params, true_stim, behaviour):
    """计算负对数似然"""
    p_common = params[0]
    n_samples = len(true_stim)
    log_likelihood = 0
    
    for i in range(n_samples):
        # 计算模型预测
        prior = compute_prior_array(x, p_common)
        likelihood = my_gaussian(x, true_stim[i], 1)
        posterior = likelihood * prior
        posterior /= np.sum(posterior)
        model_pred = np.sum(x * posterior)
        
        # 计算对数似然
        log_likelihood += -0.5 * (behaviour[i] - model_pred)**2
    
    return -log_likelihood

# 执行参数估计
true_stim, behaviour = generate_data(100, 0.95)
result = minimize(neg_log_likelihood, [0.5], args=(true_stim, behaviour), bounds=[(0, 1)])
estimated_p_common = result.x[0]

结果分析

通过上述步骤，我们可以：

1. 可视化模型组件（先验、似然、后验）的关系
2. 生成模拟参与者的行为数据
3. 从行为数据中恢复原始模型参数

这种方法可以推广到更复杂的贝叶斯模型，是认知建模中的核心技能。

总结

本教程详细介绍了：

1. 贝叶斯认知模型的构建方法
2. 如何将模型组件实现为可计算的数学形式
3. 模型反演的基本流程
4. 从行为数据中估计模型参数的完整过程

掌握这些技能对于理解人类认知过程和心理物理学实验分析至关重要。