Simulating ODEs and SDEs

0. 本 lab 导读

Lab 1 是整门课里代码量最小、但每一行的"信息密度"最高的一次作业。五个 TODO 加起来实际要写的有效代码大约 12 行，可是这 12 行覆盖了第一讲全部的核心对象：向量场（vector field）、常微分方程（ODE）、布朗运动（Brownian motion）、随机微分方程（SDE）、得分函数（score）、朗之万动力学（Langevin dynamics），以及把它们从纸上搬到计算机里的两个数值格式——欧拉法（Euler method）和欧拉–丸山法（Euler–Maruyama method）。

正因为代码短，这个 lab 的难点完全不在"写"，而在"知道自己写对了没有"。ODE 那一行错了，你会立刻看出轨迹跑飞；但 SDE 那一行如果把 $\sqrt{h}$ 写成 $h$，程序照样跑、图照样画得出来、轨迹照样有毛刺，只是物理是错的——而且步长取得越细，错得越离谱。所以本页除了逐题给参考实现，还专门用一整节讲如何用闭式统计量给自己的实现做体检：布朗运动的方差应该等于 $\sigma^2 t$，末端分布的峰度（kurtosis）应该等于 3，OU 过程的平稳方差应该等于 $\frac{\sigma^2}{2\theta}$。这些数都是能笔算出来的，跑一遍就知道有没有问题。

1. Part 0：先把接口和张量形状约定看懂

做题之前先花五分钟读懂 notebook 给的抽象基类，因为本 lab 一多半的调试时间都消耗在形状不匹配上，而这些形状约定全写在 docstring 里。

notebook 定义了两组抽象类。第一组描述"方程本身"：

class ODE(ABC):
    @abstractmethod
    def drift_coefficient(self, xt: torch.Tensor, t: torch.Tensor) -> torch.Tensor:
        # xt: (bs, dim)   t: ()  标量
        # return: (bs, dim)
        pass

class SDE(ABC):
    @abstractmethod
    def drift_coefficient(self, xt, t): ...      # (bs, dim) -> (bs, dim)
    @abstractmethod
    def diffusion_coefficient(self, xt, t): ...  # (bs, dim) -> (bs, dim)

第二组描述"怎么把方程跑起来"：

class Simulator(ABC):
    @abstractmethod
    def step(self, xt, t, dt): ...     # 你要填的就是这个

    @torch.no_grad()
    def simulate(self, x, ts):         # ts: (nts,)
        for t_idx in range(len(ts) - 1):
            t = ts[t_idx]
            h = ts[t_idx + 1] - ts[t_idx]   # 步长由 ts 相邻差自动给出
            x = self.step(x, t, h)
        return x                       # (bs, dim)

有四个约定值得提前记住，它们直接决定你 step 该怎么写：

2. Part 1：两个数值积分格式

这一部分对应 notebook 的 Question 1.1，要填两个 step。它们是全 lab 的地基：后面 Part 2、Part 3 的所有图，都是拿这两个函数跑出来的。如果这里写错，后面所有实验的结论都会是错的，而且错得很隐蔽。

Q1.1a：EulerSimulator.step

题目在问什么

给定 ODE $\frac{\dd{}}{\dd{t}}X_t = \vft{t}(X_t)$，已知当前状态 $X_t$、当前时刻 $t$、步长 $h$，求一步之后的近似状态 $X_{t+h}$。用一句数学话说：实现显式欧拉格式

数学依据

这个式子不需要背，它就是一阶泰勒展开（Taylor expansion）截断。设 $X$ 关于 $t$ 二阶可微，在 $t$ 处展开：

把 ODE 本身 $\frac{\dd{}}{\dd{t}}X_t = \vft{t}(X_t)$ 代进第一项，再把 $O(h^2)$ 那一项直接扔掉，剩下的就是欧拉格式。扔掉的那一项就是局部截断误差（local truncation error），量级 $O(h^2)$。

参考实现

class EulerSimulator(Simulator):
    def __init__(self, ode: ODE):
        self.ode = ode

    def step(self, xt: torch.Tensor, t: torch.Tensor, h: torch.Tensor):
        # 显式 Euler: x_{t+h} = x_t + h * u_t(x_t)
        # xt: (bs, dim), t: (), h: ()
        return xt + self.ode.drift_coefficient(xt, t) * h   # (bs, dim)

为什么这样写

• 为什么是 self.ode.drift_coefficient(xt, t) 而不是 (xt, t + h)？ 因为这是显式欧拉：速度取在区间左端点。取右端点得到的是隐式欧拉（backward Euler），那需要解一个关于 $X_{t+h}$ 的方程，没法一行写完；取中点得到的是中点法（midpoint method），是二阶的。课程后面提到的 Heun 法则是"先用欧拉预测，再用两端速度的平均值校正"，也是二阶。
• 形状怎么对上：drift_coefficient 返回 (bs, dim)，h 是 ()，PyTorch 把标量广播到每个元素，结果仍是 (bs, dim)，与 xt 同形状，可以直接相加。
• 为什么不用原地操作（xt += ...）：simulate_with_trajectory 会把每一步的 x.clone() 存进列表。如果你在 step 里原地改 xt，理论上 clone 已经隔离了，但一旦某处忘了 clone，整条轨迹会全部变成最后一帧的值——画出来是一堆水平直线。返回新张量最省心。
• 为什么不必自己写 torch.no_grad()：基类的 simulate / simulate_with_trajectory 已经加了装饰器。但注意 Q3.1 的 LangevinSDE 内部要用 jacrev 求导，那是 torch.func 的函数式变换，不受 no_grad 影响，所以能正常工作。

Q1.1b：EulerMaruyamaSimulator.step —— 那个 √h 从哪来

题目在问什么

给定 SDE $\dd{X_t} = \vft{t}(X_t)\dd{t} + \sigma_t \dd{W_t}$，实现它的欧拉–丸山离散：

代码上只比 Q1.1a 多一项，但那一项里的 $\sqrt{h}$ 是整个 lab、乃至整门课的分水岭。如果你只想认真读懂本页的一处，就读这里。

数学依据：布朗增量的标准差按 √h 缩放

回忆 Lecture 1 里布朗运动 $(W_t)_{t\ge 0}$ 的定义，它由三条性质刻画：

1. $W_0 = 0$；
2. 增量服从高斯：对任意 $0 \le s < t$，有 $W_t - W_s \sim \N\!\left(0,\,(t-s)I_d\right)$；
3. 增量独立：不重叠时间区间上的增量互相独立。

SDE 的形式解是一个随机积分：

欧拉–丸山做的事，就是把两个积分都用左端点的常数近似：

关键在第二行的括号。由定义中的性质 2，

也就是说它的方差是 $h$，标准差是 $\sqrt{h}$。而 $\N(0, hI_d)$ 分布的样本可以写成 $\sqrt{h}\,z$，其中 $z\sim\N(0,I_d)$——这正是代码里那个 torch.sqrt(h) * torch.randn_like(xt)。

参考实现

class EulerMaruyamaSimulator(Simulator):
    def __init__(self, sde: SDE):
        self.sde = sde

    def step(self, xt: torch.Tensor, t: torch.Tensor, h: torch.Tensor):
        # x_{t+h} = x_t + h * u_t(x_t) + sigma_t * sqrt(h) * z,  z ~ N(0, I)
        drift     = self.sde.drift_coefficient(xt, t)      # (bs, dim)
        diffusion = self.sde.diffusion_coefficient(xt, t)  # (bs, dim)
        return xt + drift * h + diffusion * torch.sqrt(h) * torch.randn_like(xt)
        #      (bs,dim) + (bs,dim)*() + (bs,dim)*()*(bs,dim)  ->  (bs, dim)

为什么这样写

• torch.randn_like(xt) 而不是 torch.randn(bs, dim)：randn_like 会自动继承 xt 的 shape、dtype 和 device。notebook 里所有张量都 .to(device) 放在 GPU 上，用 torch.randn(...) 会在 CPU 上生成，然后触发 RuntimeError: Expected all tensors to be on the same device。这是本 lab 报错率最高的一行。
• 每步必须重新采样噪声。布朗增量的独立性是定义的一部分。如果把 z 提到 __init__ 里生成一次然后每步复用，得到的就不是布朗运动，而是一条被随机斜率决定的直线族。
• torch.sqrt(h) 而不是 math.sqrt(h)：h 是 0 维 tensor，math.sqrt 对它其实也能工作（会隐式调 __float__），但如果 h 在 GPU 上会触发一次同步拷贝；更重要的是万一将来 h 变成 batched 张量，math.sqrt 就直接崩了。用 torch.sqrt 或 h ** 0.5 都行。
• 噪声乘的是 diffusion 不是 self.sde.sigma：基类接口里扩散系数是 diffusion_coefficient(xt, t) 的返回值，它可以依赖于 $x$ 和 $t$。本 lab 三个 SDE 的扩散系数恰好都是常数，但 Lab 2、Lab 3 会用到随时间变化的 $\sigma_t$。绕开接口直接读属性，到下一个 lab 就要重写。
• 漂移和扩散各调用一次，先算完再组合。LangevinSDE 的 drift_coefficient 里要跑 vmap(jacrev(...))，是本 lab 最贵的操作，别在一行里重复调用两次。

3. Part 2：把 SDE 的解画出来

这一部分让你亲手实现两个最经典的 SDE，然后画 500 条轨迹去"看"随机过程长什么样。它们不是玩具：布朗运动是所有扩散模型的噪声来源，OU 过程则是 DDPM 前向加噪过程（variance-preserving SDE）的连续时间原型。

Q2.1：实现 BrownianMotion

题目在问什么

实现 SDE

即"漂移为零、扩散系数为常数 $\sigma$"的那个 SDE。要填的是 drift_coefficient 与 diffusion_coefficient 两个方法。

数学依据

对照通式 $\dd{X_t} = \vft{t}(X_t)\dd{t} + \sigma_t\dd{W_t}$，逐项读系数：

由布朗运动的定义直接可以写出解析解：从 $X_0 = 0$ 出发，

所以均值恒为 0，标准差按 $\sigma\sqrt{t}$ 增长。这条闭式解会在第 6 节被用作测试基准。

参考实现

class BrownianMotion(SDE):
    def __init__(self, sigma: float):
        self.sigma = sigma

    def drift_coefficient(self, xt: torch.Tensor, t: torch.Tensor) -> torch.Tensor:
        # dX_t = sigma * dW_t —— 纯扩散，没有确定性漂移
        return torch.zeros_like(xt)                    # (bs, dim)

    def diffusion_coefficient(self, xt: torch.Tensor, t: torch.Tensor) -> torch.Tensor:
        # 常数扩散系数，但形状必须与 xt 对齐
        return self.sigma * torch.ones_like(xt)        # (bs, dim)

为什么这样写

• 为什么不能 return 0 或 return self.sigma？这是本题的全部考点。接口约定返回 (bs, dim) 的张量。在本题里返回 Python 标量其实碰巧也能跑通——因为 xt + 0 * h 和 xt + drift*h + sigma*sqrt(h)*randn_like(xt) 靠广播都成立。但这是运气：一旦某个 SDE 的扩散系数按维度不同（各向异性），或者你在别处对返回值调用 .shape、.device、torch.where，标量立刻炸掉。更现实的是 Lab 3：那里的扩散系数是 (bs, 1) 广播到 (bs, 1, 32, 32) 的图像张量，形状约定一旦养成坏习惯就会连环出错。
• zeros_like / ones_like 而不是 torch.zeros(xt.shape)：同样是为了自动继承 device 和 dtype。torch.zeros(xt.shape) 生成在 CPU 上，跟 GPU 上的 xt 一相加就报错。
• 为什么 t 完全没用上？因为这是时齐（time-homogeneous）的 SDE，系数不显含时间。签名里保留 t 是为了统一接口；Lab 2 的条件路径里 $\sigma_t$ 会真的随时间变。

Q2.1 文字题：σ 很大 / 接近 0 时轨迹长什么样

notebook 在实现之前先问了一次直觉，实现之后又让你改 sigma 再看一次。合并作答如下。

4. Q2.2：Ornstein–Uhlenbeck 过程与它的平稳分布

Q2.2：实现 OUProcess

题目在问什么

实现 SDE

与布朗运动相比只多了一项漂移 $-\theta X_t$：一个指向原点、大小正比于离原点距离的回复力（restoring force），物理上就是弹簧，统计上叫均值回复（mean reversion）。参数 $\theta > 0$ 是回复强度，$\sigma$ 仍是噪声强度。

数学依据

逐项读系数：漂移 $\vft{t}(x) = -\theta x$，扩散 $\sigma_t = \sigma$。代码层面就这么简单。但这一题的真正内容是它的解，因为后面三道文字题和 Q3.2 全靠它。

参考实现

class OUProcess(SDE):
    def __init__(self, theta: float, sigma: float):
        self.theta = theta
        self.sigma = sigma

    def drift_coefficient(self, xt: torch.Tensor, t: torch.Tensor) -> torch.Tensor:
        # 指向原点的线性回复力：theta 越大，被拉回得越快
        return -self.theta * xt                        # (bs, dim)

    def diffusion_coefficient(self, xt: torch.Tensor, t: torch.Tensor) -> torch.Tensor:
        return self.sigma * torch.ones_like(xt)        # (bs, dim)

为什么这样写

• 负号是全部内容。-self.theta * xt 里的负号保证了这是回复力：$x>0$ 时漂移向下，$x<0$ 时漂移向上。写成 +theta * xt 得到的是 $\dd{X_t}=\theta X_t\dd{t}+\sigma\dd{W_t}$，解含 $e^{+\theta t}$，轨迹指数发散——这个错误一眼就能从图上看出来。
• 形状自动正确。-self.theta * xt 里 self.theta 是 Python float，与 (bs, dim) 张量相乘广播后仍是 (bs, dim)，无需 ones_like。而扩散系数不含 xt，必须显式用 ones_like 撑出形状——两个方法写法不一致是有原因的，不是笔误。
• 不要在这里做任何"稳定性"修补。有人会想给漂移加 clamp 防止爆炸。不需要：只要 $\theta h < 2$，显式欧拉对这个线性系统就是稳定的。notebook 的参数下 $\theta h$ 最大约 $200 \times 0.00015 = 0.03$，非常安全。

Q2.2 文字题：θ 很小 / 很大时轨迹长什么样

Q2.2 文字题：解收敛吗？收敛到点还是收敛到分布？

Q2.2 文字题：σ–D 网格图能得出什么结论

notebook 最后那段代码把 $\sigma \in \{1,2,10\}$ 与 $D = \frac{\sigma^2}{2\theta} \in \{0.25,1,4\}$ 做笛卡儿积（$\theta$ 由 $\theta = \frac{\sigma^2}{2D}$ 反解），并且把每张图的时间轴按 $1/\sigma^2$ 缩放（ts = linspace(0, simulation_time / sigma**2, 1000)）。这个缩放不是为了好看，它本身就是结论的一部分。

5. Part 3：用 SDE 变换分布

前两部分看的是"一个点怎么动"。Part 3 换视角：一整个分布怎么被 SDE 搬运。这才是生成模型真正关心的事——我们要的不是把某个特定的噪声向量变成某张特定的图，而是把噪声分布 $\simple$ 变成数据分布 $\data$。

本部分先给出两个抽象概念，再让你实现 Langevin 动力学。这两个抽象值得单独说一下，因为它们贯穿后面所有 lab。

预备：Density 与 Sampleable 为什么要分成两个类

notebook 定义了两个互相独立的抽象基类：

这个拆分不是为了工程整洁，它编码了本课的一个核心事实：真实数据分布只满足第二条。我们手上有一百万张猫图，但没有任何办法算出"这张图的概率密度是多少"，更别说它的梯度。所以 Langevin 动力学虽然优美，却不能直接用在图像上——$\nabla\log p$ 是未知的。整门课后面要干的事，就是用神经网络把这个未知的 score 学出来（score matching），或者绕过它去学向量场（flow matching）。本 lab 让你先在"两者都有"的玩具分布上把机制跑通。

基类里 score 有一个默认实现，值得看懂：

def score(self, x: torch.Tensor) -> torch.Tensor:
    x = x.unsqueeze(1)                          # (bs, 1, dim)
    score = vmap(jacrev(self.log_density))(x)   # (bs, 1, 1, 1, dim)
    return score.squeeze((1, 2, 3))             # (bs, dim)

jacrev 对 log_density 求雅可比，vmap 把它按 batch 维向量化，从而一次性并行算出 batch 里每个样本各自的梯度。为什么不直接 torch.autograd.grad(log_density(x).sum(), x)？其实也可以（求和之后再求导，正好给出逐样本梯度，因为不同样本互不耦合），但 vmap + jacrev 的写法更函数式、不需要 requires_grad，也不受外层 @torch.no_grad() 影响——而 simulate 恰好带着这个装饰器。那一串 unsqueeze / squeeze 就是在对付 jacrev 输出的"输出形状 × 输入形状"嵌套维度。这段代码是给好的，不用改，但如果你自定义 Density 子类，记得 log_density 必须返回 (bs, 1)，否则维度对不上。

Q3.1：实现 LangevinSDE

题目在问什么

给定一个能算密度的分布 $p$，实现（过阻尼）朗之万动力学

也就是：漂移 = 半个 $\sigma^2$ 乘以 score，扩散 = $\sigma$。

数学依据：为什么偏偏是 ½σ²

系数 $\frac12\sigma^2$ 不是凑出来的，它是让 $p$ 成为平稳分布的唯一选择。用 Lecture 01 的 Fokker–Planck 方程可以一步验出来。

参考实现

class LangevinSDE(SDE):
    def __init__(self, sigma: float, density: Density):
        self.sigma = sigma
        self.density = density

    def drift_coefficient(self, xt: torch.Tensor, t: torch.Tensor) -> torch.Tensor:
        # u_t(x) = (1/2) * sigma^2 * grad log p(x)
        return 0.5 * self.sigma ** 2 * self.density.score(xt)   # (bs, dim)

    def diffusion_coefficient(self, xt: torch.Tensor, t: torch.Tensor) -> torch.Tensor:
        return self.sigma * torch.ones_like(xt)                 # (bs, dim)

为什么这样写

• 调 self.density.score(xt) 而不是自己求导。基类已经用 vmap(jacrev(...)) 实现好了，返回形状就是 (bs, dim)，与 xt 一致，可以直接线性组合。自己写 autograd.grad 要处理 requires_grad 和 no_grad 上下文，纯属自找麻烦。
• self.sigma ** 2 而不是 self.sigma。漂移里是 $\sigma^2$、扩散里是 $\sigma$，两处不同。这是本题唯一的"陷阱"，也是最值得对着上面的推导确认一遍的地方。
• 为什么漂移不含 t？目标分布 $p$ 是固定的，所以 Langevin 是时齐的 SDE。这跟后面 lab 里"沿一条概率路径 $p_t$ 逐步演化"的 SDE 不同——那时 score 会变成 $\nabla\log p_t$，显含时间。
• 性能提示：score 里的 vmap(jacrev(...)) 是整个 lab 最慢的操作。graph_dynamics 默认 1000 步 × 1000 个样本，在 CPU 上要跑好几分钟。把 device 设成 cuda，或者把 num_samples 降到 500 先调通。

Q3.1 文字题：改 σ、步数、源分布、目标密度会怎样

6. Q3.2：证明 OU 过程就是高斯的 Langevin 动力学

这是全 lab 唯一一道纯数学题，没有代码，但它是把 Part 2 和 Part 3 缝在一起的那一针。题目分两问。

Q3.2 第一问：求 N(0, σ²/2θ) 的 score

题目在问什么

证明当 $p(x) = \N\!\left(0,\ \frac{\sigma^2}{2\theta}\right)$ 时，

完整推导

Q3.2 第二问：由此说明 OU 过程就是该 p 的 Langevin 动力学

题目在问什么

把第一问的 score 代进 Langevin 动力学的定义，验证得到的 SDE 与 OU 过程逐项相同。

完整推导

7. 如何验证自己的实现是对的

本 lab 的所有产出都是图，而图有个致命问题：随机过程的图"看起来对"的容忍度极高。轨迹有毛刺、末端有钟形直方图、OU 会收敛——这些现象在 $\sqrt{h}$ 写错、$\sigma$ 少平方、噪声采两次的情况下照样出现。所以别靠眼睛，靠数。

好消息是本 lab 涉及的每个对象都有闭式统计量，可以拿蒙特卡洛估计去对拍。下面这套自查只需要几十行代码，全部在 CPU 上几秒钟跑完。我把它整理成了 labs/tests/test_labs.py，Lab 1 部分共 13 项，全部通过；下表列出每一项的理论值与实测值。

七行代码测出 Euler 是不是一阶

选一个有解析解的 ODE：$\frac{\dd{x}}{\dd{t}} = -x$，$x_0 = 1$，则 $x(1) = e^{-1}$。

class Decay(ODE):
    def drift_coefficient(self, xt, t):
        return -xt                       # (bs, dim)

sim, errs = EulerSimulator(Decay()), []
for n in [100, 200, 400]:
    ts = torch.linspace(0.0, 1.0, n + 1)
    xT = sim.simulate(torch.ones(1, 1), ts)
    errs.append(abs(xT.item() - math.exp(-1.0)))
print(errs[0] / errs[1], errs[1] / errs[2])   # 期望 ~2.0, 实测 2.004 / 2.001

两个比值都应该接近 2。如果接近 4，说明你不小心实现了二阶格式（比如中点法）；如果不趋于任何常数，说明步长没有正确进入公式。这条测试之所以有力，是因为它检验的是收敛阶——一个不依赖于常数因子、无法靠"碰巧调对参数"蒙混过关的性质。

三行代码抓出 √h 与 doubly-stochastic 两个大坑

sigma, T = 1.5, 2.0
sim = EulerMaruyamaSimulator(BrownianMotion(sigma))
xT  = sim.simulate(torch.zeros(40000, 1), torch.linspace(0.0, T, 201))

print(xT.var().item(), sigma ** 2 * T)        # 4.4871 vs 4.5000
z = (xT - xT.mean()) / xT.std()
print((z ** 4).mean().item())                 # 峰度：3.016（高斯应为 3）

Langevin 的自查：用一维高斯当目标

二维五峰混合太难收敛，不适合做单元测试。换成一维高斯 $p = \N(\mu, s^2)$，它的 score 有闭式 $-\frac{x-\mu}{s^2}$，而且单峰无势垒，$T=40$ 足够混合：

class Gauss1D(Density):
    def __init__(self, mu, s): self.mu, self.s = mu, s
    def log_density(self, x):
        return (-0.5 * ((x - self.mu) / self.s) ** 2).sum(dim=-1, keepdim=True)
    def score(self, x):                       # 覆盖基类，避免 vmap/jacrev 的形状纠缠
        return -(x - self.mu) / self.s ** 2

sim = EulerMaruyamaSimulator(LangevinSDE(sigma=1.0, density=Gauss1D(2.0, 1.3)))
xT  = sim.simulate(torch.zeros(40000, 1), torch.linspace(0.0, 40.0, 8001))
print(xT.mean().item(), xT.var().item())      # 2.0110 / 1.6840，目标 2.0 / 1.6900

注意所有粒子都从 0 出发，而目标均值是 2——如果 score 的符号或系数错了，均值绝对到不了 2.0110。方差 1.6840 对 $\frac12\sigma^2$ 这个系数尤其敏感：把 0.5 * sigma ** 2 误写成 0.5 * sigma（$\sigma=1$ 时看不出来）或 sigma ** 2（漏 $\frac12$，方差会掉到目标的一半左右），这一项立刻报警。所以做这类自查时永远别用 $\sigma=1$ 和 $\mu=0$。

8. 本 lab 小结

延伸阅读

• Lecture 01 · Flow and Diffusion Models —— 本 lab 全部公式的理论出处：向量场与流、布朗运动的定义、Euler / Euler–Maruyama 的推导、连续性方程与 Fokker–Planck 方程的完整证明。
• Generative Modeling by Estimating Gradients of the Data Distribution (Song & Ermon, 2019) —— 把本 lab 的 Langevin 动力学推向真实数据：正是这篇论文发现"单一 $p$ 上的 Langevin 混合太慢"，并提出用多尺度噪声退火（annealed Langevin）来解决，直接对应本页 Q3.1 观察到的模式偏权重问题。
• Score-Based Generative Modeling through Stochastic Differential Equations (Song et al., 2021) —— 把 DDPM 与 score matching 统一到 SDE 框架；其中的 VP-SDE 就是系数随时间变化的 OU 过程，Q3.2 那道题是理解它的钥匙。
• Flow Matching for Generative Modeling (Lipman et al., 2022) —— Lab 2 的理论来源，用 ODE 而非 SDE 完成同一件事。
• Kloeden & Platen, Numerical Solution of Stochastic Differential Equations —— SDE 数值方法的标准参考书。想知道 Euler–Maruyama 的强收敛阶为何是 $\frac12$ 而弱收敛阶是 $1$（本 lab 关心的其实是后者，因为我们只在乎分布），以及 Milstein 格式如何把强阶提到 1，看这本。