格兰杰因果——神经层面正确，经过区域不同的血流动力学延迟后可能反向

# 模拟两个脑区的“神经”信号：x 驱动 y（y_t 依赖 x_{t-1}），反向耦合为 0。
# granger_F 用最小二乘拟合 p 阶向量自回归，返回 Geweke 形式的 F_{x→y} = ln(受限残差方差 / 完整残差方差)。
# 然后把两路信号分别与峰值时间不同的伽马型血流动力学响应函数（HRF）卷积并降采样，模拟 BOLD：
# 情形 1：两区 HRF 相同；情形 2：被驱动的 y 区 HRF 峰值反而更早（区域间 HRF 差异约 1–2 秒在文献中常见）。
# 输出解读：神经信号上 F_{x→y} 远大于 F_{y→x}；HRF 相同时方向大体保留；
# y 区 HRF 更早时，BOLD 上的“因果方向”可能被削弱甚至反转。数值为模拟结果，仅说明机制。
import numpy as np
from scipy.stats import gamma

def granger_F(src, dst, p=2):
    """Geweke 形式的格兰杰因果量 F_{src→dst} = ln(var_restricted / var_full)。"""
    T = len(dst)
    Y = dst[p:]
    lags_dst = np.column_stack([dst[p - k:T - k] for k in range(1, p + 1)])
    lags_src = np.column_stack([src[p - k:T - k] for k in range(1, p + 1)])
    ones = np.ones((T - p, 1))
    Xr = np.hstack([ones, lags_dst])
    Xf = np.hstack([ones, lags_dst, lags_src])
    rr = Y - Xr @ np.linalg.lstsq(Xr, Y, rcond=None)[0]
    rf = Y - Xf @ np.linalg.lstsq(Xf, Y, rcond=None)[0]
    return np.log(rr.var() / rf.var())

def hrf(t, peak):
    """简化的单伽马 HRF，峰值约在 peak 秒（形状参数 a = peak + 1，尺度 1）。"""
    h = gamma.pdf(t, a=peak + 1, scale=1.0)
    return h / h.sum()

rng = np.random.default_rng(0)
dt, T = 0.1, 60000                       # 神经信号 10 Hz、共 6000 秒（模拟数据）
x = np.zeros(T); y = np.zeros(T)
for t in range(1, T):                    # 一阶自回归神经动力学，x → y，延迟 0.1 s
    x[t] = 0.9 * x[t - 1] + rng.standard_normal()
    y[t] = 0.9 * y[t - 1] + 0.5 * x[t - 1] + rng.standard_normal()

print(f"神经信号（10 Hz）：F_x→y = {granger_F(x, y):.3f}，F_y→x = {granger_F(y, x):.3f}")

th = np.arange(0, 30, dt)
step = int(2.0 / dt)                     # 降采样到 TR = 2 s
Fres = [(granger_F(x, y), granger_F(y, x))]   # 记录各情形的 F 值，供作图
for name, px, py in [("两区 HRF 峰值相同(5 s, 5 s)", 5.0, 5.0),
                     ("y 区 HRF 更早(x: 6 s, y: 4.5 s)", 6.0, 4.5)]:
    bx = np.convolve(x, hrf(th, px))[:T][::step]
    by = np.convolve(y, hrf(th, py))[:T][::step]
    bx = bx + 0.05 * bx.std() * rng.standard_normal(len(bx))   # 少量测量噪声
    by = by + 0.05 * by.std() * rng.standard_normal(len(by))
    fxy, fyx = granger_F(bx, by), granger_F(by, bx)
    verdict = "x→y" if fxy > fyx else "y→x"
    print(f"BOLD（TR=2 s），{name}：F_x→y = {fxy:.3f}，F_y→x = {fyx:.3f}，判定方向 {verdict}")
    Fres.append((fxy, fyx))

import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
for pk, ls, lab in [(5.0, "-", "peak 5 s (equal case)"), (6.0, "--", "x region, peak 6 s"), (4.5, ":", "y region, peak 4.5 s")]:
    axes[0].plot(th, hrf(th, pk) / dt, ls, lw=2, label=lab)
axes[0].set(xlabel="Time (s)", ylabel="HRF (a.u.)", title="Regional HRFs used in the simulation")
axes[0].legend(fontsize=8)
Fres = np.array(Fres)
pos = np.arange(3)
axes[1].bar(pos - 0.2, Fres[:, 0], 0.4, color="#4477AA", label="F x->y (true direction)")
axes[1].bar(pos + 0.2, Fres[:, 1], 0.4, color="#CC6677", label="F y->x")
axes[1].set_xticks(pos, ["Neural (10 Hz)", "BOLD, equal HRF", "BOLD, y HRF earlier"], fontsize=9)
axes[1].set(ylabel="Geweke Granger F", title="HRF differences can flip the inferred direction")
axes[1].legend(fontsize=8)
plt.tight_layout(); plt.show()