Diffusionモデルの拡散過程の確率微分方程式を
d \boldsymbol{X}_t = f(t) \boldsymbol{X}_t dt + g(t) d \boldsymbol{W}_t
とする。$\boldsymbol{X}_t$は$d$次元の確率過程で、$\boldsymbol{W}_t$はWiener過程。
通常$f(t) \le 0$で
f(t) \boldsymbol{X}_t dt
の項は$\boldsymbol{X}_t$が平均的に$\boldsymbol{0}$に近づいていくことを表し、
g(t) d \boldsymbol{W}_t
の項は$t$から$t+dt$の間で加えられるノイズが分散$g(t)^2dt$程度の各次元が独立な正規分布であることを表す。
時刻$0$における$\boldsymbol{X}_0$の分布（初期分布）は既知であるとして、各時刻$t$の$\boldsymbol{X}_t$分布を求める。
まず、ノイズの項のない常微分方程式
d \boldsymbol{x} = f(t) \boldsymbol{x} dt
は変数分離型の一階の常微分方程式なので、簡単に解くことができる。$\boldsymbol{x}$の時刻$0$での初期値を$\boldsymbol{x}_0$とすれば、
\boldsymbol{x} = e^{\int_0^{t} f(\xi) d\xi} \boldsymbol{x}_0
である。
次に$d$次元の確率過程$\boldsymbol{Y}_t$を
\boldsymbol{Y}_t = e^{-\int_0^{t} f(\xi) d\xi} \boldsymbol{X}_t
と定義する。時刻$0$においては、$\boldsymbol{Y}_0 = \boldsymbol{X}_0$である。
$h(t, \boldsymbol{x}) = e^{-\int_0^{t} f(\xi) d\xi} \boldsymbol{x}$とすると、$\boldsymbol{Y}_t = h(t, \boldsymbol{X}_t)$であり、
\begin{eqnarray}
\frac{\partial h}{\partial t} &=& -f(t) e^{-\int_0^{t} f(\xi) d\xi} \boldsymbol{x} \\
\frac{\partial h}{\partial x_j} &=& e^{-\int_0^{t} f(\xi) d\xi} \\
\frac{\partial^2 h}{\partial x_i \partial x_j} &=& 0 \\
\end{eqnarray}
なので、伊藤の公式から、
d \boldsymbol{Y}_t = e^{-\int_0^{t} f(\xi) d\xi} g(t) d \boldsymbol{W}_t
が成り立つ。
この確率過程の時刻$t>0$の分布は積分から求められる。
\begin{eqnarray}
\boldsymbol{Y}_t - \boldsymbol{Y}_0 &=& \int_0^t e^{-\int_0^{\tau} f(\xi) d\xi} g(\tau) d \boldsymbol{W}_{\tau} \\
&=& \lim_{\Delta \rightarrow 0} \sum_{j=1}^{N} e^{-\int_0^{t_j} f(\xi) d\xi} g(t_j) (\boldsymbol{W}_{t_{j+1}} - \boldsymbol{W}_{t_j})
\end{eqnarray}
ただし、$t_1=0, t_N=t, \Delta = \sup_j(t_{j+1} - t_j)$。
$\boldsymbol{W}_t$はWiener過程なので、
\boldsymbol{W}_{t_{j+1}} - \boldsymbol{W}_{t_j} \sim N(\boldsymbol{0}, (t_{j+1} - t_j)\mathbf{I}_d)
であり、$\boldsymbol{W}_2 - \boldsymbol{W}_1$, $\boldsymbol{W}_3 - \boldsymbol{W}_2$, $...$は独立である。ただし$\mathbf{I}_d$は$d$次元単位行列。
正規分布に従う確率変数の和はまた正規分布であり、分散は和となるため、
\sum_{j=1}^{N} e^{-\int_0^{t_j} f(\xi) d\xi} g(t_j) (\boldsymbol{W}_{t_{j+1}} - \boldsymbol{W}_{t_j}) \sim N \left(\boldsymbol{0}, \sum_{j=1}^{N} e^{-2\int_0^{t_j} f(\xi) d\xi} g(t_j)^2 (t_{j+1} - t_j) \mathbf{I}_d\right)
である。
\lim_{\Delta \rightarrow 0} \sum_{j=1}^{N} e^{- 2 \int_0^{t_j} f(\xi) d\xi} g(t_j)^2 (t_{j+1} - t_j) = \int_0^t e^{-\int_0^{\tau} f(\xi) ds} g(\tau) d\tau
であるから、
\boldsymbol{Y}_t - \boldsymbol{Y}_0 \sim N \left(\boldsymbol{0}, \int_0^t e^{-2\int_0^{\tau} f(\xi) d\xi} g(\tau)^2 d\tau \mathbf{I}_d \right)
となる。
\begin{eqnarray}
s(t) &=& e^{\int_0^{t} f(\xi) d\xi} \\
\sigma(t) &=& \sqrt{\int_0^t \frac{g(\tau)^2}{s(\tau)^2} d\tau}
\end{eqnarray}
とすれば、$\boldsymbol{Y}_t - \boldsymbol{Y}_0 \sim N \left(\boldsymbol{0}, \sigma (t)^2 \mathbf{I}_d\right)$である。
$\boldsymbol{Y}_0 = \boldsymbol{X}_0$なので、$\boldsymbol{X}_0$の確率密度関数を$p_0(\boldsymbol{x})$とすれば、これは$\boldsymbol{Y}_0$の確率密度関数でもある。$\boldsymbol{Y}_t$の確率密度関数を$p_Y(\boldsymbol{y}, t)$とすると、$\boldsymbol{Y}_t - \boldsymbol{Y}_0 \sim N \left(\boldsymbol{0}, \sigma (t)^2 \mathbf{I}_d\right)$と$\boldsymbol{Y}_0 \sim p_0$から、
\begin{eqnarray}
p_Y(\boldsymbol{y}, t) &=& \left[ p_0 * N(\boldsymbol{0}, \sigma (t)^2 \mathbf{I}_d)\right] (\boldsymbol{y})\\
&=& \int_{\boldsymbol{z} \in \mathbb{R}^d} p_0(\boldsymbol{z}) \mathcal{N}_d \left(\boldsymbol{y} - \boldsymbol{z}\middle| \boldsymbol{0}, \sigma (t)^2 \mathbf{I}_d \right) d \boldsymbol{z}
\end{eqnarray}
である。ただし、$[\cdot * \cdot]$は畳み込み、$\mathcal{N}_d(\cdot | \mu, \Sigma)$は平均$\mu$、分散$\Sigma$の正規分布の確率密度関数を表している。
$\boldsymbol{Y}_t = \frac{1}{s(t)} \boldsymbol{X}_t$なので、$\boldsymbol{X}_t$の確率密度関数$p_X(\boldsymbol{x}, t)$とすると、
\begin{eqnarray}
p_X(\boldsymbol{x}, t) &=& \frac{1}{s(t)^d} p_Y \left( \frac{\boldsymbol{x}}{s(t)}, t \right) \\
&=& \frac{1}{s(t)^d} \int_{\boldsymbol{z} \in \mathbb{R}^d} p_0(\boldsymbol{z}) \mathcal{N}_d\left(\frac{\boldsymbol{x}}{s(t)} - \boldsymbol{z}\middle| \boldsymbol{0}, \sigma(t)^2 \mathbf{I}_d \right) d \boldsymbol{z} \\
&=& \int_{\boldsymbol{z} \in \mathbb{R}^d} p_0(\boldsymbol{z}) \mathcal{N}_d\left(\boldsymbol{x} - s(t) \boldsymbol{z}\middle| \boldsymbol{0}, s(t)^2 \sigma(t)^2 \mathbf{I}_d \right) d \boldsymbol{z} \\
&=& \int_{\boldsymbol{z} \in \mathbb{R}^d} p_0(\boldsymbol{z}) \mathcal{N}_d\left(\boldsymbol{x}\middle| s(t) \boldsymbol{z}, s(t)^2 \sigma(t)^2 \mathbf{I}_d \right) d \boldsymbol{z} \\
\end{eqnarray}
となる。これが確率微分方程式$d \boldsymbol{X}_t = f(t) \boldsymbol{X}_t dt + g(t) d \boldsymbol{W}_t$の時刻$0$での初期分布の確率密度関数が$p_0(\boldsymbol{x})$の場合の時刻$t$における分布の確率密度関数である。
この確率密度関数の式から、時刻$0$において$\boldsymbol{x}_0$にある場合、時刻$t$における$\boldsymbol{X}_t$の条件付き確率分布は正規分布
N \left(s(t) \boldsymbol{x}_0, s(t)^2 \sigma(t)^2 \mathbf{I}_d \right)
であることがわかる。
また、導出の途中式の
p_X(\boldsymbol{x}, t) = \frac{1}{s(t)^d} \int_{\boldsymbol{z} \in \mathbb{R}^d} p_0(\boldsymbol{z}) \mathcal{N}_d\left(\frac{\boldsymbol{x}}{s(t)} - \boldsymbol{z}\middle| \boldsymbol{0}, \sigma(t)^2 \mathbf{I}_d \right) d \boldsymbol{z}
より、畳み込みで
p_X(\boldsymbol{x}, t) = \frac{1}{s(t)^d}\left[ p_{\mathrm{data}} * N(\boldsymbol{0}, \sigma(t)^2 \mathbf{I}_d)\right] \left( \frac{\boldsymbol{x}}{s(t)} \right)
と表現できる。$\boldsymbol{X}_t$の分布は、初期分布$p_0$に分散$\sigma(t)^2$の大きさのノイズを加えて、スケールを$\frac{1}{s(t)}$にしたものである。