重積分の変数変換とヤコビ行列式 — なぜ|det J|が現れるのか

円板上の積分を計算するとき、私たちはほとんど反射的に $dx\,dy \to r\,dr\,d\theta$ と書き換えます。球なら $dx\,dy\,dz \to r^2\sin\theta\,dr\,d\theta\,d\varphi$ です。ではこの $r$ や $r^2\sin\theta$ は、いったいどこから湧いて出てきたのでしょうか。「公式だから」と暗記して済ませることもできますが、この係数の正体を知らないまま先に進むと、少し違う座標変換を持ち出された瞬間に手が止まります。

この係数の正体は ヤコビ行列式（ヤコビアン）$\det \bm{J}$ の絶対値 であり、その意味はたった一言で言えます — 「変換が局所的に体積を何倍に引き伸ばしているか」 です。円板を極座標で測るとき、半径 $r$ の位置では $\theta$ 方向に $1$ 進むと実際には弧長 $r$ だけ進むので、面積が $r$ 倍に膨らむ。それだけの話です。

左はもとの座標のマス目で、どれも同じ大きさです。これを変換すると（右）、マス目は場所ごとに広がったり縮んだりします。色は各マス目の面積が何倍になったかを実測したもので、いちばん広がった場所は $1.13$ 倍、いちばん縮んだ場所は $0.64$ 倍と、同じ図の中で倍率が2倍近く違っています。この「場所ごとの面積の倍率」を数式で表したものが、これから見るヤコビ行列式です。

この「体積の伸縮率」という視点を手に入れると、一見バラバラだった話が一本につながります。

• ガウス積分 $\int_{-\infty}^{\infty}e^{-x^2}dx=\sqrt{\pi}$ — 原始関数が初等関数で書けないこの積分が、極座標に移った瞬間に暗算レベルで解けてしまう理由がわかります。正規分布の正規化定数 $1/\sqrt{2\pi}$ の出所でもあります。
• 確率密度の変換 — 確率変数を $\bm{Y}=g(\bm{X})$ と変換したとき密度がどう変わるか。多変量正規分布の密度に $\sqrt{\det\bm{\Sigma}}$ が現れる理由も、実はヤコビアンです。
• 正規化フロー（Normalizing Flow） — 深層生成モデルの一族は、まさにこの変数変換公式を対数の形で書き下したものを損失関数にしています。「ヤコビ行列式が速く計算できるニューラルネットを設計する」という、一見奇妙な研究テーマの動機がここにあります。

本記事の内容

• 線形写像が体積を $|\det \bm{A}|$ 倍することを、平行四辺形・平行六面体の体積公式から導く
• 微分可能写像は局所的に線形であるという事実から、変数変換公式をリーマン和で組み立てる
• 極座標の $r$、球座標の $r^2\sin\theta$ を1行ずつ計算し、ガウス積分を極座標で解く
• 確率密度の変換公式 $p_Y(\bm{y})=p_X(g^{-1}(\bm{y}))\,|\det \bm{J}_{g^{-1}}(\bm{y})|$ と正規化フローへの接続
• Pythonで「面積比 = $|\det\bm{J}|$」をメッシュで実測し、モンテカルロ積分と密度変形で公式を数値検証する

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

• ヤコビ行列とは？定義と意味をわかりやすく解説
• 行列式とは？定義と幾何学的な意味
• 重積分とは？累次積分の考え方
• 球座標系（極座標系）の基礎

すでに公式の使い方だけを知りたい方は ヤコビアンと変数変換（重積分の座標変換公式） に各座標系の計算例がまとまっています。本記事はその「なぜ」の部分を掘り下げる立場で書いています。

1変数の置換積分を思い出すと、疑問が浮かび上がる

高校で習った置換積分を書いてみます。$x=g(t)$ と置くと

$$ \int_{g(a)}^{g(b)} f(x)\,dx = \int_a^b f(g(t))\,g'(t)\,dt $$

でした。右辺に $g'(t)$ という余計な因子が付きます。学校では「$dx = g'(t)\,dt$ だから」と説明されますが、これを幾何的に読み直すと次の意味になります。$t$ 軸上の微小区間 $dt$ は、写像 $g$ によって $x$ 軸上の長さ $g'(t)\,dt$ の区間に引き伸ばされる。 つまり $g'(t)$ は「局所的な長さの伸縮率」です。積分は「値 × 幅」の足し算ですから、幅が $g’$ 倍になるなら、その分を掛けておかないと総和が合いません。

下側の横軸に置いた2つの区間は、どちらも同じ幅 $\Delta t = 0.28$ です。ところが左の縦軸に写された像の幅は、傾きが緩い場所では $0.29$、傾きが急な場所では $0.59$ と2倍以上違います。その比 $0.29/0.28 \approx 1.04$、$0.59/0.28 \approx 2.12$ が、まさにその点での $g'(t)$ の値になっています。積分するときにこの倍率を掛け忘れれば、右側の領域の寄与を半分に見積もってしまうことになります。

ここで自然な疑問が生まれます。2変数、3変数ではどうなるのか。$(u,v)$ 平面の微小な正方形が $(x,y)$ 平面に写されるとき、その「面積の伸縮率」は何で書けるのでしょうか。1変数のときは伸縮率がスカラー $g’$ ひとつで済みましたが、2変数では写像は $u$ 方向にも $v$ 方向にも、しかも斜めにも伸び縮みできます。偏微分は $\partial x/\partial u,\ \partial x/\partial v,\ \partial y/\partial u,\ \partial y/\partial v$ の4つあり、これらを「面積の倍率」というひとつの数にまとめる必要があります。

その「まとめ方」が行列式です。以下では、まず線形写像という一番簡単な場合で「行列式 = 体積の倍率」を確認し、その後に一般の曲がった写像へ話を広げます。

線形写像は体積を $|\det \bm{A}|$ 倍する

最初に相手にするのは、いちばん素直な写像 — 線形写像です。

$$ \bm{x} = \bm{A}\bm{u}, \qquad \bm{A} \in \mathbb{R}^{n\times n} $$

イメージを掴むために、$2\times 2$ の場合で単位正方形がどこへ行くかを考えます。$(u,v)$ 平面の単位正方形の頂点は $(0,0), (1,0), (1,1), (0,1)$ です。これらを $\bm{A}$ で写すと、それぞれ $\bm{0}$、第1列 $\bm{a}_1$、$\bm{a}_1+\bm{a}_2$、第2列 $\bm{a}_2$ に移ります。つまり 単位正方形は、$\bm{A}$ の列ベクトルが張る平行四辺形に変形します。 元の面積は $1$ でしたから、平行四辺形の面積そのものが「面積の倍率」になります。

2次元：平行四辺形の面積が $|\det\bm{A}|$ になること

$\bm{a}_1=(a_{11}, a_{21})^\top$、$\bm{a}_2=(a_{12},a_{22})^\top$ とします。2辺の長さと挟む角 $\theta$ を使えば、平行四辺形の面積は

$$ S = \|\bm{a}_1\|\,\|\bm{a}_2\|\,\sin\theta $$

です。$\sin\theta$ が扱いにくいので、$\sin^2\theta = 1-\cos^2\theta$ と内積 $\bm{a}_1\cdot\bm{a}_2 = \|\bm{a}_1\|\|\bm{a}_2\|\cos\theta$ を使って書き換えます。両辺を2乗すると

$$ S^2 = \|\bm{a}_1\|^2\|\bm{a}_2\|^2\left(1-\cos^2\theta\right) = \|\bm{a}_1\|^2\|\bm{a}_2\|^2 – (\bm{a}_1\cdot\bm{a}_2)^2 $$

となり、三角関数が消えて成分だけの式になりました。ここに成分を代入して展開します。

$$ S^2 = (a_{11}^2+a_{21}^2)(a_{12}^2+a_{22}^2) – (a_{11}a_{12}+a_{21}a_{22})^2 $$

まず前半を展開すると $a_{11}^2a_{12}^2 + a_{11}^2a_{22}^2 + a_{21}^2a_{12}^2 + a_{21}^2a_{22}^2$、後半を展開すると $a_{11}^2a_{12}^2 + 2a_{11}a_{12}a_{21}a_{22} + a_{21}^2a_{22}^2$ です。引き算すると $a_{11}^2a_{12}^2$ と $a_{21}^2a_{22}^2$ がきれいに相殺して

$$ S^2 = a_{11}^2a_{22}^2 – 2a_{11}a_{22}a_{12}a_{21} + a_{12}^2a_{21}^2 = (a_{11}a_{22}-a_{12}a_{21})^2 $$

が残ります。右辺は完全平方の形です。$S \ge 0$ なので平方根をとって

$$ \begin{equation} S = |a_{11}a_{22}-a_{12}a_{21}| = |\det \bm{A}| \end{equation} $$

を得ます。「$2\times2$ 行列式 $ad-bc$」という、意味のわからない暗記事項だったものが、単位正方形が写った先の面積そのもの だったわけです。

3枚のパネルはいずれも、単位正方形（灰色）が $\bm{A}$ の2本の列ベクトルで張られる平行四辺形に写る様子です。中央では靴ひも公式で測った実測面積 $2.50$ が $\det\bm{A}=2.50$ とぴたりと一致しています。右は $x$ と $y$ を入れ替える行列で、面積は $1.00$ のまま変わらないのに $\det\bm{A}=-1.00$ と符号だけが負になっています。頂点番号 $1\to2\to3\to4$ の回り方が反時計回りから時計回りへ反転していることが、この符号の正体です。

3次元：平行六面体の体積とスカラー三重積

3次元では、単位立方体は $\bm{A}$ の3本の列ベクトル $\bm{a}_1,\bm{a}_2,\bm{a}_3$ が張る平行六面体に写ります。平行六面体の体積は「底面積 × 高さ」で

$$ V = \underbrace{\|\bm{a}_1\times\bm{a}_2\|}_{\text{底面積}} \times \underbrace{\left|\bm{a}_3\cdot \frac{\bm{a}_1\times\bm{a}_2}{\|\bm{a}_1\times\bm{a}_2\|}\right|}_{\text{高さ（底面法線方向の成分）}} = \left|(\bm{a}_1\times\bm{a}_2)\cdot\bm{a}_3\right| $$

と書けます。ここで外積 $\bm{a}_1\times\bm{a}_2$ の大きさが底面（$\bm{a}_1,\bm{a}_2$ が張る平行四辺形）の面積であること、そして高さは $\bm{a}_3$ の「底面の単位法線方向の成分」であることを使いました。最後に現れた $(\bm{a}_1\times\bm{a}_2)\cdot\bm{a}_3$ は スカラー三重積 と呼ばれ、成分で書くとちょうど $3\times3$ 行列式のサラスの展開になります。

$$ (\bm{a}_1\times\bm{a}_2)\cdot\bm{a}_3 = a_{13}(a_{21}a_{32}-a_{31}a_{22}) – a_{23}(a_{11}a_{32}-a_{31}a_{12}) + a_{33}(a_{11}a_{22}-a_{21}a_{12}) = \det\bm{A} $$

したがって $V = |\det \bm{A}|$ です。2次元と同じ結論が3次元でも成り立ちました。

3本のベクトル $\bm{a}_1,\bm{a}_2,\bm{a}_3$ が作る箱が、単位立方体の行き先です。塗った面が底面（面積 $2.37$）、破線が $\bm{a}_3$ から底面へ下ろした高さ（$1.28$）で、その積 $3.03$ が $|\det\bm{A}|=3.03$ と一致しています。$\bm{a}_3$ を底面と平行な方向にいくら動かしても高さは変わらない — これは行列式が「ある列に他の列の定数倍を足しても不変」という性質そのものを、絵で見ていることになります。

$n$ 次元：「体積らしさ」を要求すると行列式しかない

$n$ 次元では外積が使えませんが、逆に「体積とはどういう性質を持つべきか」から攻めると、行列式が必然的に出てきます。$n$ 本のベクトルに対して符号付き体積を返す関数 $D(\bm{a}_1,\dots,\bm{a}_n)$ に、次の3つを要求します。

1. 各引数について線形（多重線形性） — 1辺を2倍に伸ばせば体積も2倍。1辺を2つのベクトルの和に分ければ体積も分配される。
2. 2つの引数が等しいと $0$（交代性） — 辺が重なれば平行体は潰れて体積 $0$。
3. $D(\bm{e}_1,\dots,\bm{e}_n)=1$（正規化） — 単位立方体の体積は $1$。

線形代数の基本定理として、この3条件を満たす関数は行列式ただひとつです。つまり 「体積」という概念を素直に公理化すると、行列式が一意に定まる。だから行列式が体積倍率として現れるのは偶然ではありません。

符号の意味と、なぜ絶対値をつけるのか

$\det\bm{A}$ は負にもなります。負のときは何が起きているかというと、向き（orientation）が反転 しています。2次元なら鏡映（裏返し）、3次元なら右手系が左手系になる状態です。たとえば $\bm{A}=\begin{pmatrix}0&1\\1&0\end{pmatrix}$（$x$ と $y$ の入れ替え）は $\det\bm{A}=-1$ で、面積は変わらず向きだけが反転します。

積分で使う「面積・体積」は非負の量（測度）なので、向きの情報は捨てて $|\det\bm{A}|$ を使います。これが変数変換公式に絶対値が付く理由です。逆に、微分形式を使って向き付き積分を扱う立場では絶対値を外し、符号を積分の向きに押し付けます。1変数の置換積分で絶対値が見当たらないのは、まさにこの流儀 — $g$ が減少関数のときは $g(a)>g(b)$ となって積分区間の向きが逆転し、その向きが符号を吸収しているからです。

ここまでは写像が線形、つまり「まっすぐ」な場合の話でした。しかし極座標変換のような現実の変換は曲がっています。次はこのギャップを埋めます。

微分可能写像は局所的に線形 — だからヤコビ行列式が効く

極座標変換 $(r,\theta)\mapsto(r\cos\theta, r\sin\theta)$ は、格子を扇状にぐにゃりと曲げます。到底「線形写像」とは言えません。しかし 虫めがねで拡大していくと、曲がった写像もどんどん線形写像に見えてきます。 地球は丸いのに、目の前の地面は平らに見えるのと同じ理屈です。

これを数式にしたのが、多変数版のテイラー展開（全微分可能性）です。写像 $g:\mathbb{R}^n\to\mathbb{R}^n$ が点 $\bm{u}_0$ で微分可能なら

$$ \begin{equation} g(\bm{u}_0 + \Delta\bm{u}) = g(\bm{u}_0) + \bm{J}_g(\bm{u}_0)\,\Delta\bm{u} + o(\|\Delta\bm{u}\|) \end{equation} $$

が成り立ちます。ここで $\bm{J}_g$ は ヤコビ行列

$$ \bm{J}_g(\bm{u}) = \begin{pmatrix} \dfrac{\partial x_1}{\partial u_1} & \cdots & \dfrac{\partial x_1}{\partial u_n}\\[6pt] \vdots & \ddots & \vdots\\[6pt] \dfrac{\partial x_n}{\partial u_1} & \cdots & \dfrac{\partial x_n}{\partial u_n} \end{pmatrix} $$

です。式(2)が言っているのは、「$\bm{u}_0$ の近くでは、$g$ は平行移動 $g(\bm{u}_0)$ と線形写像 $\bm{J}_g(\bm{u}_0)$ の合成にすぎない」 ということ。誤差項 $o(\|\Delta\bm{u}\|)$ は $\Delta\bm{u}$ より高位の微小量なので、十分小さい領域では無視できます。

平行移動は体積を変えません。線形写像 $\bm{J}_g(\bm{u}_0)$ は前節の結果から体積を $|\det\bm{J}_g(\bm{u}_0)|$ 倍します。したがって:

$\bm{u}_0$ のまわりの微小な立方体（体積 $h^n$）は、$g$ によって体積がおよそ $|\det\bm{J}_g(\bm{u}_0)|\,h^n$ の微小平行体に写る。

これが「ヤコビ行列式 = 局所的な体積伸縮率」の中身です。$\det\bm{J}$ は ヤコビアン とも呼ばれ、$\dfrac{\partial(x_1,\dots,x_n)}{\partial(u_1,\dots,u_n)}$ と書かれることもあります。

同じ点のまわりを、窓の半幅を $0.35 \to 0.12 \to 0.03$ と狭めながら見た図です。緑の実線が実際の像、黒の破線がヤコビ行列による線形近似（平行四辺形）です。左では両者がはっきりずれていて格子も湾曲していますが、右ではほとんど重なって区別がつきません。面積比も $0.9310 \to 0.9201 \to 0.9186$ と、$\det\bm{J}=0.9185$ へ着実に近づいています。「拡大すれば線形」という主張が、そのまま目に見えるかたちで現れています。

リーマン和で変数変換公式を組み立てる

局所の話が済んだので、これを足し合わせて積分にします。$g: U \to X$ を全単射な $C^1$ 級写像とし、$U$ を一辺 $h$ の小立方体 $\{Q_k\}$ に分割します。各小立方体の代表点を $\bm{u}_k$ とすると、その像 $g(Q_k)$ は

$$ \mathrm{vol}\big(g(Q_k)\big) \approx |\det \bm{J}_g(\bm{u}_k)|\,h^n $$

の微小平行体です。いま $X$ 上の積分 $\int_X f(\bm{x})\,d\bm{x}$ をリーマン和で近似すると、小片 $g(Q_k)$ 上で $f$ はほぼ一定値 $f(g(\bm{u}_k))$ をとるので

$$ \int_X f(\bm{x})\,d\bm{x} \approx \sum_k f\big(g(\bm{u}_k)\big)\,\mathrm{vol}\big(g(Q_k)\big) $$

となります。ここに先ほどの体積の近似式を代入すると

$$ \int_X f(\bm{x})\,d\bm{x} \approx \sum_k f\big(g(\bm{u}_k)\big)\,|\det \bm{J}_g(\bm{u}_k)|\,h^n $$

が得られます。この右辺をよく見てください。これは $U$ 上の関数 $f(g(\bm{u}))\,|\det\bm{J}_g(\bm{u})|$ のリーマン和そのもの です。$h\to 0$ の極限をとれば、両辺はそれぞれ本物の積分に収束して

$$ \begin{equation} \int_{X} f(\bm{x})\,d\bm{x} = \int_{U} f\big(g(\bm{u})\big)\,\left|\det \bm{J}_g(\bm{u})\right|\,d\bm{u} \end{equation} $$

が成り立ちます。これが 重積分の変数変換公式 です。証明としてはラフで、「誤差項の和が本当に $0$ に潰れるか」を詰めるのが厳密な証明の山場になります（誤差は各小片で $o(h^n)$ ですが、小片の個数は $O(h^{-n})$ 個あるため、単純な足し算では消えてくれません。実際には $C^1$ 級の仮定から誤差が一様に効くことを使います）。ただし公式が持つ意味 — 「被積分関数に局所体積伸縮率を掛けて、変換前の座標で積分し直す」 — はこの議論で完全に見えています。

定理の正確なステートメント

使うときに引っかかりやすい条件を明示しておきます。

定理（変数変換） $U, X \subset \mathbb{R}^n$ を開集合、$g: U \to X$ を $C^1$ 級の全単射で、すべての $\bm{u}\in U$ で $\det\bm{J}_g(\bm{u})\neq 0$ とする。このとき可積分な $f$ に対して式(3)が成り立つ。

実用上の注意は3つです。

1. 全単射性は「測度ゼロの例外」を許してよい。 極座標変換は $r=0$ で単射になりませんし（原点に全ての $\theta$ が写る）、$\theta=0$ と $\theta=2\pi$ も同じ点です。しかしこれらは面積 $0$ の集合なので、積分値に影響しません。だから安心して $r\in[0,\infty),\ \theta\in[0,2\pi)$ を使えます。
2. $\det\bm{J}=0$ の点でも同様。 極座標の原点では $\det\bm{J}=r=0$ ですが、これも1点なので問題ありません。逆に、$\det\bm{J}$ が広い領域でゼロなら、その変換は次元を潰しており変数変換として使えません。
3. どちら向きの写像か混同しない。 式(3)の $g$ は「新しい変数 $\bm{u}$ から元の変数 $\bm{x}$ への写像」です。極座標なら $g(r,\theta)=(r\cos\theta, r\sin\theta)$ で、$(r,\theta)$ が新変数です。向きを間違えると $|\det\bm{J}|$ が逆数になります。

逆写像のヤコビアンは逆数になる

向きの話が出たので、逆写像との関係を押さえておきます。$g$ が可逆なら $g^{-1}\circ g = \mathrm{id}$ で、連鎖律を適用すると

$$ \bm{J}_{g^{-1}}\big(g(\bm{u})\big)\,\bm{J}_g(\bm{u}) = \bm{I} $$

つまり $\bm{J}_{g^{-1}}(g(\bm{u})) = \bm{J}_g(\bm{u})^{-1}$ です。両辺の行列式をとり、積の行列式は行列式の積という性質 $\det(\bm{A}\bm{B})=\det\bm{A}\det\bm{B}$ を使うと

$$ \begin{equation} \det \bm{J}_{g^{-1}}\big(g(\bm{u})\big) = \frac{1}{\det \bm{J}_g(\bm{u})} \end{equation} $$

を得ます。片方向に体積が $c$ 倍になるなら、戻る向きでは $1/c$ 倍になる という、当たり前すぎるほど当たり前の関係です。この式は後で確率密度の変換を扱うときにそのまま効いてきます。

理屈は揃いました。ここからは実際に手を動かして、極座標と球座標のヤコビアンを計算してみます。

極座標のヤコビアン — $r$ の正体

極座標変換は

$$ x = r\cos\theta,\qquad y = r\sin\theta $$

でした。新変数は $(r,\theta)$、旧変数は $(x,y)$ です。ヤコビ行列を書くには4つの偏微分を並べます。

$$ \frac{\partial x}{\partial r} = \cos\theta,\quad \frac{\partial x}{\partial \theta} = -r\sin\theta,\quad \frac{\partial y}{\partial r} = \sin\theta,\quad \frac{\partial y}{\partial \theta} = r\cos\theta $$

これを並べると

$$ \bm{J} = \begin{pmatrix} \cos\theta & -r\sin\theta \\ \sin\theta & r\cos\theta\end{pmatrix} $$

$2\times2$ の行列式 $ad-bc$ を計算します。

$$ \det \bm{J} = \cos\theta \cdot r\cos\theta – (-r\sin\theta)\cdot\sin\theta = r\cos^2\theta + r\sin^2\theta $$

共通因子 $r$ でくくり、$\cos^2\theta+\sin^2\theta=1$ を使うと

$$ \begin{equation} \det \bm{J} = r\left(\cos^2\theta+\sin^2\theta\right) = r \end{equation} $$

よって $dx\,dy = r\,dr\,d\theta$ です。

幾何で確かめる

いまの計算を、微小な扇形の面積で検算してみましょう。半径 $r$ から $r+dr$、角度 $\theta$ から $\theta+d\theta$ で囲まれる領域は、2つの扇形の差です。半径 $R$、中心角 $d\theta$ の扇形の面積は $\frac{1}{2}R^2 d\theta$ ですから

$$ dS = \frac{1}{2}(r+dr)^2 d\theta – \frac{1}{2}r^2 d\theta = \frac{1}{2}\left(2r\,dr + (dr)^2\right)d\theta $$

$(dr)^2$ は $dr$ に比べて高位の微小量なので落とすと

$$ dS \approx r\,dr\,d\theta $$

行列式の計算と完全に一致しました。式(5)の $r$ は、「半径が大きいほど、同じ角度幅でも弧が長くなる」 という素朴な事実を表しているだけなのです。

左の $(r,\theta)$ 平面ではすべてのセルが同じ大きさの長方形ですが、中央の $(x,y)$ 平面に写すと外側ほど大きな扇形になります。同じ $\Delta r = 0.30$、$\Delta\theta=0.35$ でも、内側のセルの面積は $0.047$、外側は $0.205$ と4倍以上の差がつきました。右のグラフは各セルの面積を実測して半径に対してプロットしたもので、原点を通る直線にきれいに乗っています。この直線の傾きが $\Delta r\,\Delta\theta$、つまり比例係数として現れているのが $|\det\bm{J}|=r$ です。$r=0$（原点）では $\det\bm{J}=0$ となり、$\theta$ をどう動かしても同じ点に写ってしまう — 変換が潰れることも数式に正しく反映されています。

この $r$ が1個あるおかげで、次に見る有名な積分が解けてしまいます。

ガウス積分を極座標で解く

$$ I = \int_{-\infty}^{\infty} e^{-x^2}\,dx $$

この積分の被積分関数 $e^{-x^2}$ には初等関数で書ける原始関数が存在しません。まともに攻めると詰みます。ところが、「2乗して2次元にする」 という一手で状況が一変します。

まず $I$ を2つ掛け合わせます。積分変数はダミーなので、片方を $y$ と書き直せます。

$$ I^2 = \left(\int_{-\infty}^{\infty} e^{-x^2}dx\right)\left(\int_{-\infty}^{\infty} e^{-y^2}dy\right) $$

それぞれの積分は互いに独立なので、フビニの定理により1つの重積分にまとめられます。

$$ I^2 = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} e^{-x^2}e^{-y^2}\,dx\,dy = \iint_{\mathbb{R}^2} e^{-(x^2+y^2)}\,dx\,dy $$

ここで指数の肩に現れた $x^2+y^2$ に注目してください。これは 原点からの距離の2乗、つまり被積分関数は回転対称です。回転対称な関数を扱うなら、座標も回転対称なものにするのが自然 — 極座標の出番です。式(3)に $g(r,\theta)=(r\cos\theta,r\sin\theta)$、$|\det\bm{J}|=r$ を代入すると

$$ I^2 = \int_0^{2\pi}\!\!\int_0^{\infty} e^{-r^2}\;\underbrace{r\,dr\,d\theta}_{=\,dx\,dy} $$

被積分関数は $\theta$ を含まないので、$\theta$ 積分は単に $2\pi$ を出すだけです。

$$ I^2 = 2\pi \int_0^{\infty} r\,e^{-r^2}\,dr $$

残った積分に $s=r^2$ と置換します。$ds = 2r\,dr$ すなわち $r\,dr = \frac{1}{2}ds$ で、積分区間は $s: 0 \to \infty$ です。

$$ I^2 = 2\pi \int_0^{\infty} \frac{1}{2}e^{-s}\,ds = \pi\left[-e^{-s}\right]_0^{\infty} = \pi\,(0-(-1)) = \pi $$

$I>0$（正値関数の積分）なので

$$ \begin{equation} I = \int_{-\infty}^{\infty} e^{-x^2}\,dx = \sqrt{\pi} \approx 1.7724539 \end{equation} $$

が得られました。

この解法の核心はどこにあったのでしょうか。 決め手は、ヤコビアンとして現れた $r$ が、$\int r\,e^{-r^2}dr$ の形をつくったことです。もし $r$ が無ければ $\int e^{-r^2}dr$ のままで、1次元のときと同じく手が出ません。$r$ という因子が付いたおかげで、$e^{-r^2}$ の微分がちょうど $-2r\,e^{-r^2}$ であることと噛み合い、置換一発で解けたのです。「変数変換で余分な因子が付いて面倒」ではなく、「余分な因子こそが問題を解いてくれる」 — これがヤコビアンの面白いところです。

左は $e^{-(x^2+y^2)}$ の曲面で、等高線が真円になる — つまり原点からの距離だけで値が決まることが見て取れます。右は動径方向の被積分関数の比較で、青が $e^{-r^2}$、赤がヤコビアンを掛けた $r\,e^{-r^2}$ です。青は $r=0$ で最大値 $1$ を取りますが、赤は $0$ から立ち上がって $r=1/\sqrt{2}$ で山を作ります。この赤い曲線の下の面積が厳密に $1/2$ であり、$\theta$ 積分の $2\pi$ を掛けて $I^2=\pi$ が出ます。原点で $0$ になるのは「原点近傍は面積が小さいので寄与しない」という当たり前の事実の反映です。

なお、$\sigma$ 付きの一般形は $x \to x/(\sqrt{2}\sigma)$ の置換で $\int_{-\infty}^{\infty}e^{-x^2/(2\sigma^2)}dx = \sqrt{2\pi}\sigma$ となり、正規分布の密度関数の正規化定数 $1/(\sqrt{2\pi}\sigma)$ が出てきます。詳しくは ガウス積分の導出をわかりやすく解説 を参照してください。

2次元の様子がわかったところで、3次元に上げてみましょう。計算量は増えますが、やることは同じです。

球座標のヤコビアン — $r^2\sin\theta$ を余因子展開で導く

球座標（物理系の慣習で、$\theta$ を天頂角、$\varphi$ を方位角とします）は

$$ x = r\sin\theta\cos\varphi,\qquad y = r\sin\theta\sin\varphi,\qquad z = r\cos\theta $$

で定義されます。$r\in[0,\infty),\ \theta\in[0,\pi],\ \varphi\in[0,2\pi)$ です。9個の偏微分を計算します。

$r$ で微分すると（$r$ は1次で入っているので、単に $r$ を落とすだけ）:

$$ \frac{\partial x}{\partial r}=\sin\theta\cos\varphi,\qquad \frac{\partial y}{\partial r}=\sin\theta\sin\varphi,\qquad \frac{\partial z}{\partial r}=\cos\theta $$

$\theta$ で微分すると（$\sin\theta \to \cos\theta$、$\cos\theta \to -\sin\theta$）:

$$ \frac{\partial x}{\partial \theta}=r\cos\theta\cos\varphi,\qquad \frac{\partial y}{\partial \theta}=r\cos\theta\sin\varphi,\qquad \frac{\partial z}{\partial \theta}=-r\sin\theta $$

$\varphi$ で微分すると（$z$ は $\varphi$ を含まないのでゼロ）:

$$ \frac{\partial x}{\partial \varphi}=-r\sin\theta\sin\varphi,\qquad \frac{\partial y}{\partial \varphi}=r\sin\theta\cos\varphi,\qquad \frac{\partial z}{\partial \varphi}=0 $$

並べると

$$ \bm{J} = \begin{pmatrix} \sin\theta\cos\varphi & r\cos\theta\cos\varphi & -r\sin\theta\sin\varphi\\ \sin\theta\sin\varphi & r\cos\theta\sin\varphi & r\sin\theta\cos\varphi\\ \cos\theta & -r\sin\theta & 0 \end{pmatrix} $$

第3行に $0$ があるので、第3行に沿って余因子展開する のが得策です。第3行の成分は $(\cos\theta,\ -r\sin\theta,\ 0)$ で、余因子の符号は $(+,-,+)$ の並びです。

$$ \det\bm{J} = \cos\theta \cdot M_{31} – (-r\sin\theta)\cdot M_{32} + 0 $$

ここで $M_{31}, M_{32}$ は小行列式です。まず $M_{31}$（第3行・第1列を除いた $2\times2$）を計算します。

$$ M_{31} = \begin{vmatrix} r\cos\theta\cos\varphi & -r\sin\theta\sin\varphi \\ r\cos\theta\sin\varphi & r\sin\theta\cos\varphi\end{vmatrix} = r^2\sin\theta\cos\theta\cos^2\varphi + r^2\sin\theta\cos\theta\sin^2\varphi $$

$r^2\sin\theta\cos\theta$ でくくって $\cos^2\varphi+\sin^2\varphi=1$ を使うと

$$ M_{31} = r^2\sin\theta\cos\theta $$

次に $M_{32}$（第3行・第2列を除いた $2\times2$）です。

$$ M_{32} = \begin{vmatrix} \sin\theta\cos\varphi & -r\sin\theta\sin\varphi \\ \sin\theta\sin\varphi & r\sin\theta\cos\varphi\end{vmatrix} = r\sin^2\theta\cos^2\varphi + r\sin^2\theta\sin^2\varphi = r\sin^2\theta $$

これらを展開式に代入します。

$$ \det\bm{J} = \cos\theta\cdot r^2\sin\theta\cos\theta + r\sin\theta\cdot r\sin^2\theta = r^2\sin\theta\cos^2\theta + r^2\sin^3\theta $$

$r^2\sin\theta$ でくくると、またしても三角関数の基本関係が効いて

$$ \begin{equation} \det\bm{J} = r^2\sin\theta\left(\cos^2\theta+\sin^2\theta\right) = r^2\sin\theta \end{equation} $$

$\theta\in[0,\pi]$ では $\sin\theta \ge 0$ なので絶対値は不要で、体積要素は

$$ dx\,dy\,dz = r^2\sin\theta\,dr\,d\theta\,d\varphi $$

となります。

$r^2\sin\theta$ の幾何的な読み方

この式は2つの因子の積として読めます。

• $r^2$ — 半径方向に離れるほど、同じ立体角に対応する球面上の面積が $r^2$ で増える。太陽から遠ざかると光が弱まる「逆2乗則」と同じ理屈です。
• $\sin\theta$ — 方位角 $\varphi$ 方向の弧長は、天頂角 $\theta$ の位置では半径 $r\sin\theta$ の円周上を進むので、$\sin\theta$ の重みが付く。地球儀の緯線を思い浮かべてください。赤道（$\theta=\pi/2$、$\sin\theta=1$）の緯線は一番長く、極（$\theta=0,\pi$、$\sin\theta=0$）では点に潰れます。

検算として半径 $a$ の球の体積を計算しましょう。

$$ V = \int_0^{2\pi}\!\!\int_0^{\pi}\!\!\int_0^{a} r^2\sin\theta\,dr\,d\theta\,d\varphi $$

3つの積分が完全に分離しているので順に計算できます。$\varphi$ 積分は $2\pi$、$\theta$ 積分は $\int_0^\pi \sin\theta\,d\theta = [-\cos\theta]_0^\pi = 2$、$r$ 積分は $\int_0^a r^2 dr = a^3/3$ です。

$$ V = 2\pi \cdot 2 \cdot \frac{a^3}{3} = \frac{4}{3}\pi a^3 $$

見慣れた公式が出ました。球の体積公式の $4/3$ は、$\theta$ 積分の $2$ と $r$ 積分の $1/3$ の産物 だったわけです。

左の球面では、$\theta$ と $\varphi$ を等間隔に刻んだマス目が赤道付近で広く、極付近では細長く潰れています。中央のグラフはその理由で、単位半径での緯線の長さ $2\pi\sin\theta$ は赤道で最大 $2\pi$、両極でちょうど $0$ になります。右は $r^2\sin\theta$ を中点則でリーマン和にした結果で、分割数 $N$ を増やすと誤差が $N^{-2}$ で落ちながら $4\pi/3 = 4.18879$ に収束しています。体積要素の式が正しいことの数値的な裏付けです。

座標変換への応用はここまでです。次は視点を変えて、この同じ公式が確率論でどう使われるかを見ます。

確率密度の変数変換 — 密度は「体積あたり」の量である

確率変数 $\bm{X}$ が密度 $p_X$ を持つとき、これを $\bm{Y}=g(\bm{X})$ と変換したら $\bm{Y}$ の密度はどうなるでしょうか。「$p_X$ に $g^{-1}$ を代入すればいい」と思いたくなりますが、それでは足りません。

直感から入りましょう。確率密度は「確率質量 ÷ 体積」です。 ある小領域の中にある確率の量（質量）は、変換しても変わりません — 元の領域に入っていた粒子は、写った先の領域にそのまま入っているからです。しかし 入れ物である領域の体積は $|\det\bm{J}|$ 倍に変わります。 質量が同じで体積が2倍になれば、密度は半分。だから密度には $1/|\det\bm{J}_g|$ という因子、すなわち $|\det\bm{J}_{g^{-1}}|$ が掛かるのです。人口が同じ市が合併で面積2倍になれば人口密度が半分になる、という話とまったく同じ構造です。

導出

$g$ を可逆な $C^1$ 級写像とします。任意の可測集合 $B$ に対し、$\bm{Y}\in B$ と $\bm{X}\in g^{-1}(B)$ は同じ事象なので

$$ P(\bm{Y}\in B) = P\big(\bm{X}\in g^{-1}(B)\big) = \int_{g^{-1}(B)} p_X(\bm{x})\,d\bm{x} $$

右辺に変数変換公式(3)を適用します。今回は「新変数 $\bm{y}$ から旧変数 $\bm{x}$ への写像」が $g^{-1}$ なので、$\bm{x}=g^{-1}(\bm{y})$ と置いて $d\bm{x} = |\det\bm{J}_{g^{-1}}(\bm{y})|\,d\bm{y}$ とします。積分領域 $g^{-1}(B)$ は $\bm{y}$ の世界では $B$ に戻ります。

$$ P(\bm{Y}\in B) = \int_{B} p_X\big(g^{-1}(\bm{y})\big)\,\left|\det \bm{J}_{g^{-1}}(\bm{y})\right|\,d\bm{y} $$

一方、定義から $P(\bm{Y}\in B) = \int_B p_Y(\bm{y})\,d\bm{y}$ です。この等式が 任意の $B$ で成り立つのですから、被積分関数どうしが（ほとんど至るところで）一致していなければなりません。よって

$$ \begin{equation} p_Y(\bm{y}) = p_X\big(g^{-1}(\bm{y})\big)\,\left|\det \bm{J}_{g^{-1}}(\bm{y})\right| \end{equation} $$

を得ます。式(4)の逆数関係を使えば、$\bm{x}=g^{-1}(\bm{y})$ として

$$ \begin{equation} p_Y(\bm{y}) = \frac{p_X(\bm{x})}{\left|\det \bm{J}_{g}(\bm{x})\right|} \end{equation} $$

とも書けます。この形のほうが「体積が $|\det\bm{J}_g|$ 倍に膨らんだ分だけ密度が薄まる」という直感に素直に対応しています。実装では両辺の対数をとった

$$ \begin{equation} \log p_Y(\bm{y}) = \log p_X(\bm{x}) – \log\left|\det \bm{J}_{g}(\bm{x})\right| \end{equation} $$

が使われます。掛け算が引き算になり、数値的にも安定します。

例1：多変量正規分布の $\sqrt{\det\bm{\Sigma}}$ はヤコビアンだった

$\bm{Z}\sim\mathcal{N}(\bm{0},\bm{I})$、つまり密度 $p_Z(\bm{z}) = (2\pi)^{-n/2}\exp(-\|\bm{z}\|^2/2)$ から出発して、アフィン変換 $\bm{Y} = \bm{A}\bm{Z}+\bm{\mu}$（$\bm{A}$ は正則）を考えます。ヤコビ行列は定数行列 $\bm{A}$ そのものなので $\det\bm{J}_g = \det\bm{A}$、逆写像は $\bm{z}=\bm{A}^{-1}(\bm{y}-\bm{\mu})$ です。式(9)に代入すると

$$ p_Y(\bm{y}) = \frac{1}{(2\pi)^{n/2}|\det\bm{A}|}\exp\left(-\frac{1}{2}\left\|\bm{A}^{-1}(\bm{y}-\bm{\mu})\right\|^2\right) $$

指数の中身を整理します。ノルムの2乗は内積なので

$$ \left\|\bm{A}^{-1}(\bm{y}-\bm{\mu})\right\|^2 = (\bm{y}-\bm{\mu})^\top (\bm{A}^{-1})^\top\bm{A}^{-1}(\bm{y}-\bm{\mu}) = (\bm{y}-\bm{\mu})^\top (\bm{A}\bm{A}^\top)^{-1}(\bm{y}-\bm{\mu}) $$

ここで共分散行列 $\bm{\Sigma}=\bm{A}\bm{A}^\top$ とおきます。さらに $\det\bm{\Sigma} = \det\bm{A}\cdot\det\bm{A}^\top = (\det\bm{A})^2$ なので $|\det\bm{A}| = \sqrt{\det\bm{\Sigma}}$ です。代入すると

$$ p_Y(\bm{y}) = \frac{1}{(2\pi)^{n/2}\sqrt{\det\bm{\Sigma}}}\exp\left(-\frac{1}{2}(\bm{y}-\bm{\mu})^\top\bm{\Sigma}^{-1}(\bm{y}-\bm{\mu})\right) $$

多変量正規分布の密度そのものです。分母の $\sqrt{\det\bm{\Sigma}}$ は「標準正規分布を $\bm{A}$ で引き伸ばしたときの体積膨張率」に他ならない ことがわかりました。1次元で $1/\sigma$ が付くのの自然な一般化です。

例2：2次元正規分布を極座標に移すとレイリー分布が出る

$(X,Y)$ が独立に $\mathcal{N}(0,1)$ に従うとき、その同時密度は $p(x,y)=\frac{1}{2\pi}e^{-(x^2+y^2)/2}$ です。これを極座標 $(R,\Theta)$ に移します。$|\det\bm{J}|=r$ を掛けて

$$ p(r,\theta) = \frac{1}{2\pi}e^{-r^2/2}\cdot r $$

$\theta$ について $0$ から $2\pi$ まで積分すれば（被積分関数は $\theta$ に依存しないので $2\pi$ 倍するだけ）、$R$ の周辺密度は

$$ p_R(r) = r\,e^{-r^2/2},\qquad r\ge 0 $$

これは レイリー分布 です。原点からの距離の分布が原点で $0$ になる（$p_R(0)=0$）のは奇妙に見えるかもしれませんが、これこそヤコビアン $r$ の効果です。原点近傍は面積が小さいので、密度が高くても確率質量が集まらない。$r$ を掛けることで、この「面積の重み」が正しく反映されているわけです。ちなみにレイリー分布は無線通信のマルチパスフェージングのモデルとしても現れ、その導出はまさにこの計算そのものです。

左の散布図で、原点に近い半径 $0.25$〜$0.60$ の輪（面積 $0.93$）に入る標本は全体の $13.4\%$、外側の半径 $1.50$〜$1.85$ の輪（面積 $3.68$）に入るのは $14.4\%$ です。外側の輪は密度が数分の1しかないのに、面積が4倍あるおかげで同じくらいの確率を集めています。右のグラフでは、$r$ を掛け忘れた $e^{-r^2/2}$（破線）が原点で最大になるのに対し、正しいレイリー分布（実線）は原点で $0$ から立ち上がり $r=1$ で山を作ります。ヒストグラムは実線とぴたり重なっており、面積の重み $r$ が確かに効いていることがわかります。

例3：対数正規分布

$X\sim\mathcal{N}(0,1)$ に対し $Y=e^X$ とします。1次元なのでヤコビアンはスカラー $dy/dx = e^x = y$ です。逆写像は $x=\log y$、$|dx/dy| = 1/y$ なので式(8)より

$$ p_Y(y) = \frac{1}{\sqrt{2\pi}}\exp\left(-\frac{(\log y)^2}{2}\right)\cdot\frac{1}{y},\qquad y>0 $$

これが対数正規分布です。後ほどPythonで数値検証します。

密度の変換がわかると、その応用として現代の生成モデルが視野に入ってきます。

正規化フロー — 変数変換公式を損失関数にする

深層生成モデルで「尤度を厳密に計算したい」という要求は根強くあります。VAEは変分下界（近似）しか出せず、GANはそもそも尤度を定義しません。これに対し 正規化フロー（Normalizing Flow） は、変数変換公式をそのまま使って厳密な対数尤度を計算します。

考え方は単純です。単純な基底分布（普通は標準正規分布）$p_Z$ からのサンプル $\bm{z}$ を、可逆なニューラルネット $\bm{x} = G_\theta(\bm{z})$ で複雑な分布に変形する。逆向きに $f_\theta = G_\theta^{-1}$ とおけば、データ $\bm{x}$ の対数尤度は式(10)から

$$ \begin{equation} \log p_\theta(\bm{x}) = \log p_Z\big(f_\theta(\bm{x})\big) + \log\left|\det \bm{J}_{f_\theta}(\bm{x})\right| \end{equation} $$

と書けます。右辺は完全に計算可能な量なので、これを最大化するように $\theta$ を学習すればよい — 近似も下界もありません。

ただし、素朴に実装すると壁にぶつかります。$n$ 次元データに対し $\bm{J}_{f_\theta}$ は $n\times n$ 行列で、その行列式の計算は一般に $O(n^3)$ です。画像なら $n$ は数千〜数万ですから、1サンプルごとにこれを計算するのは非現実的です。

そこで正規化フローの研究は、「表現力を保ちつつ、ヤコビ行列式が安く計算できるネットワークをどう設計するか」 という一点に集中してきました。代表的なトリックが アフィンカップリング層 です。入力を2つに分け、片方はそのまま通し、もう片方を片方の関数でスケール・シフトします。

$$ \bm{y}_{1:d} = \bm{x}_{1:d},\qquad \bm{y}_{d+1:n} = \bm{x}_{d+1:n}\odot\exp\big(s(\bm{x}_{1:d})\big) + t(\bm{x}_{1:d}) $$

このときヤコビ行列は

$$ \bm{J} = \begin{pmatrix} \bm{I}_d & \bm{0}\\ \ast & \mathrm{diag}\big(\exp(s(\bm{x}_{1:d}))\big)\end{pmatrix} $$

という 下三角ブロック行列 になります。三角行列の行列式は対角成分の積なので、左下の $\ast$ がどれだけ複雑でも一切計算する必要がありません。

$$ \log|\det\bm{J}| = \sum_{i=1}^{n-d} s_i(\bm{x}_{1:d}) $$

$O(n^3)$ が $O(n)$ の足し算になりました。しかも $s,t$ は任意に複雑なニューラルネットでよく、逆変換も $\bm{x}_{d+1:n} = (\bm{y}_{d+1:n}-t)\odot\exp(-s)$ と解析的に書けます。RealNVPやGlowといった手法はこの層を積み重ねて構成されます。

「行列式の形が計算コストを決め、それがネットワークアーキテクチャを決める」— 重積分の変数変換という古典的な定理が、現代の深層生成モデルの設計原理を直接規定しているのは、なかなか痛快な話ではないでしょうか。詳細は 正規化フロー（Normalizing Flow）とは と RealNVPのアフィンカップリング層 を参照してください。

理論は一通り揃いました。最後に、ここまでの主張を数値で確かめます。

Pythonでの実装と検証

3つの実験を行います。(1) 微小セルの面積比が本当に $|\det\bm{J}|$ になるかをメッシュで実測する、(2) モンテカルロ積分で変数変換公式が成り立つことを確認する（ついでに $|\det\bm{J}|$ を忘れるとどう間違うかを見る）、(3) 可逆変換で確率密度がどう変形するかを可視化する、です。

実験1：単位正方形の像と面積比

まず、はっきり曲がった写像を用意します。

$$ g(u,v) = \big(u + 0.3\sin 2v,\ \ v + 0.3\sin 2u\big) $$

ヤコビ行列は

$$ \bm{J}_g = \begin{pmatrix} 1 & 0.6\cos 2v\\ 0.6\cos 2u & 1\end{pmatrix},\qquad \det\bm{J}_g = 1 – 0.36\cos 2u\cos 2v $$

です。$|\cos| \le 1$ より $\det\bm{J}_g \in [0.64, 1.36]$ で、単位正方形上で常に正 — つまり折り返しの起きない可逆な変形です。この写像で単位正方形の格子を写してみます。

import numpy as np
import matplotlib, matplotlib.pyplot as plt

# 日本語フォント設定
for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
    if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
        plt.rcParams["font.family"] = cand
        break
plt.rcParams["axes.unicode_minus"] = False

def g(u, v):
    """(u,v) -> (x,y) の非線形な可逆写像"""
    return u + 0.3*np.sin(2*v), v + 0.3*np.sin(2*u)

def detJ(u, v):
    """ヤコビ行列式 det J = 1 - 0.36 cos(2u) cos(2v)"""
    return 1 - 0.36*np.cos(2*u)*np.cos(2*v)

# 格子線を描くための細かいサンプル点
t = np.linspace(0, 1, 200)
lines = np.linspace(0, 1, 11)

fig, axes = plt.subplots(1, 2, figsize=(11, 5))
for c in lines:
    axes[0].plot(t, np.full_like(t, c), color="tab:blue", lw=0.9)
    axes[0].plot(np.full_like(t, c), t, color="tab:red", lw=0.9)
    xh, yh = g(t, np.full_like(t, c))          # v = c の像
    xv, yv = g(np.full_like(t, c), t)          # u = c の像
    axes[1].plot(xh, yh, color="tab:blue", lw=0.9)
    axes[1].plot(xv, yv, color="tab:red", lw=0.9)

axes[0].set_title("変換前：$(u,v)$ 平面の単位正方形と格子")
axes[0].set_xlabel("$u$"); axes[0].set_ylabel("$v$")
axes[1].set_title("変換後：$(x,y)$ 平面での像（格子が曲がる）")
axes[1].set_xlabel("$x$"); axes[1].set_ylabel("$y$")
for ax in axes:
    ax.set_aspect("equal"); ax.grid(alpha=0.3)
plt.tight_layout()
plt.show()

左のまっすぐな格子が、右では波打った網目に変形しています（各マスを面積の倍率で塗り分けたものが、冒頭に掲げた概念図です）。注目してほしいのは 網目の「マス目の大きさ」が場所によって違う ことです。左では全部のマスが同じ大きさなのに、右では膨らんでいる場所と縮んでいる場所があります。この膨らみ具合こそが $|\det\bm{J}|$ で、これから数値で確かめます。

次に、単位正方形を $N\times N$ のセルに分割し、各セルの像の面積を実測して $|\det\bm{J}|$ と比べます。像は4頂点を結んだ四角形で近似し、面積は靴ひも公式（多角形の座標から面積を求める公式）で求めます。

import numpy as np

def cell_area_ratio(N):
    """単位正方形を N×N に分割し、(セル像の面積)/(元の面積) と det J を比較"""
    e = np.linspace(0, 1, N + 1)
    h = 1.0 / N
    ratios, dets = [], []
    for i in range(N):
        for j in range(N):
            u0, u1, v0, v1 = e[i], e[i+1], e[j], e[j+1]
            # セルの4頂点の像（反時計回り）
            P = [g(u0, v0), g(u1, v0), g(u1, v1), g(u0, v1)]
            x = np.array([p[0] for p in P]); y = np.array([p[1] for p in P])
            # 靴ひも公式で四角形の面積
            A = 0.5*abs(np.dot(x, np.roll(y, -1)) - np.dot(y, np.roll(x, -1)))
            ratios.append(A / h**2)
            dets.append(detJ(u0, v0))   # セル左下隅で評価
    return np.array(ratios), np.array(dets)

print(f"{'N':>5} {'最大相対誤差':>14} {'平均相対誤差':>14}")
for N in [4, 8, 16, 32, 64, 128]:
    r, d = cell_area_ratio(N)
    rel = np.abs(r - d) / d
    print(f"{N:5d} {rel.max():14.5f} {rel.mean():14.6f}")

実行結果は次のようになります。

    N   最大相対誤差   平均相対誤差
    4        0.12245       0.079030
    8        0.05853       0.037135
   16        0.02885       0.018029
   32        0.01424       0.008869
   64        0.00708       0.004398
  128        0.00353       0.002190

この表は理論を見事に裏書きしています。$N$ を2倍にする（セルの一辺 $h$ を半分にする）たびに、誤差がきれいに半分になっています。 つまり誤差は $O(h)$ で、$h\to 0$ でゼロに収束します。これは式(2)のテイラー展開で捨てた項が $o(\|\Delta\bm{u}\|)$ だったことの直接の反映です。$N=4$ では12%もずれていた面積比が、$N=128$ では0.35%まで一致する — 「微分可能写像は局所的に線形」という主張が、数値としてこの目で見えました。

なお、セルの左下隅ではなく中心で $\det\bm{J}$ を評価すると誤差は $O(h^2)$ に改善します（$N$ を2倍にすると誤差が1/4）。中点則が台形則より高精度なのと同じ理屈で、片側評価の偏りが打ち消し合うためです。

続いて、面積の総和がリーマン和として積分に収束することも確認しておきます。

import numpy as np

exact = 1 - 0.09*np.sin(2)**2   # ∫∫ det J du dv の解析値

def riemann_sum(N):
    e = np.linspace(0, 1, N + 1)[:-1]   # 各セルの左下隅
    h = 1.0 / N
    U, V = np.meshgrid(e, e, indexing="ij")
    return np.sum(detJ(U, V)) * h * h

print(f"解析値 ∫∫|det J| du dv = {exact:.8f}")
for N in [4, 8, 16, 32, 64, 128]:
    S = riemann_sum(N)
    print(f"N={N:4d}  リーマン和={S:.8f}  誤差={abs(S-exact):.2e}")

出力は次の通りです。

解析値 ∫∫|det J| du dv = 0.92558604
N=   4  リーマン和=0.86065210  誤差=6.49e-02
N=   8  リーマン和=0.89471761  誤差=3.09e-02
N=  16  リーマン和=0.91060696  誤差=1.50e-02
N=  32  リーマン和=0.91821727  誤差=7.37e-03
N=  64  リーマン和=0.92193273  誤差=3.65e-03
N= 128  リーマン和=0.92376726  誤差=1.82e-03

こちらも誤差が $N$ の倍増ごとに半減し、解析値 $1-0.09\sin^2 2 = 0.92558604$ に収束しています。ここでの解析値は $\int_0^1\!\!\int_0^1 (1-0.36\cos 2u\cos 2v)\,du\,dv = 1 – 0.36\left(\frac{\sin 2}{2}\right)^2$ を手計算したものです。「像の総面積 = $\int\int|\det\bm{J}|$」 という式(3)の $f\equiv 1$ の場合が、数値で確認できました。

2つの表を両対数グラフにしたものです。左では、セルの隅で $\det\bm{J}$ を評価した最大・平均相対誤差がどちらも $N^{-1}$ の破線と平行に並んでおり、誤差が $O(h)$ であることがひと目でわかります。セル中心で評価した緑の系列だけは傾きが倍の $N^{-2}$ になっていて、$N=128$ では平均相対誤差が $10^{-5}$ 台まで落ちています。右のリーマン和の誤差も同じく $N^{-1}$ に乗り、解析値 $0.92558604$ に向かって単調に減っています。「$h\to0$ で誤差が消える」という極限操作が、絵として確認できたことになります。

実験2：モンテカルロ積分で変換公式を検証

次は、変数変換公式が積分値を本当に保つのかを確かめます。題材は半径 $2$ の円板 $D$ 上のガウス型関数の積分です。

$$ \iint_D e^{-(x^2+y^2)}\,dx\,dy,\qquad D=\{(x,y): x^2+y^2\le 4\} $$

極座標に移すと $\int_0^{2\pi}\!\!\int_0^2 e^{-r^2}r\,dr\,d\theta = 2\pi\cdot\frac{1}{2}(1-e^{-4}) = \pi(1-e^{-4})$ と厳密に解けます。この厳密値を基準に、(a) 直交座標での棄却法、(b) 極座標で $r$ を掛けた版、(c) 極座標で $r$ を掛け忘れた 版の3つを比較します。

import numpy as np

rng = np.random.default_rng(0)
M = 2_000_000
exact = np.pi * (1 - np.exp(-4))

# (a) 直交座標：[-2,2]^2 に一様サンプルし、円内だけ拾う
x = rng.uniform(-2, 2, M); y = rng.uniform(-2, 2, M)
inside = x**2 + y**2 <= 4
I_cart = 16.0 * np.mean(np.where(inside, np.exp(-(x**2 + y**2)), 0.0))

# (b) 極座標：(r,θ) に一様サンプルし、|det J| = r を掛ける
r  = rng.uniform(0, 2, M)
th = rng.uniform(0, 2*np.pi, M)
I_polar = (2 * 2*np.pi) * np.mean(np.exp(-r**2) * r)

# (c) 誤り例：ヤコビアン r を掛け忘れた場合
I_wrong = (2 * 2*np.pi) * np.mean(np.exp(-r**2))

print(f"厳密値                       : {exact:.6f}")
print(f"(a) 直交座標 MC              : {I_cart:.6f}  相対誤差 {abs(I_cart-exact)/exact:.2e}")
print(f"(b) 極座標 MC（r を掛ける）  : {I_polar:.6f}  相対誤差 {abs(I_polar-exact)/exact:.2e}")
print(f"(c) 極座標 MC（r を忘れる）  : {I_wrong:.6f}  ← 全く違う値になる")

実行結果は次のようになります。

厳密値                       : 3.084052
(a) 直交座標 MC              : 3.087545  相対誤差 1.13e-03
(b) 極座標 MC（r を掛ける）  : 3.084967  相対誤差 2.97e-04
(c) 極座標 MC（r を忘れる）  : 5.548034  ← 全く違う値になる

3つの数値から読み取れることが3点あります。第一に、(a) と (b) は同じ値に収束しています。 座標系がまったく違うのに同じ積分値が出る — 変数変換公式が成立している何よりの証拠です。第二に、(c) はヤコビアン $r$ を落としただけで $5.55$ という無関係な数になりました。 理論値 $\pi^{3/2}\mathrm{erf}(2)=5.542$ に収束しており、これは「別の積分を計算してしまった」ことを意味します。$r$ は飾りではなく、積分の値そのものを決めている因子なのです。第三に、この試行では (b) の誤差が (a) の約4分の1 に収まりました。直交座標では標本の $1-\pi/4\approx 21\%$ が円の外に落ちて無駄になるうえ、指示関数の不連続がばらつきを生みます。極座標なら全標本が有効に使えるので分散が小さい。1回の試行では乱数の引きに左右されるので、後で乱数の種を変えて60回繰り返した平均も確認しますが、そこでも極座標のほうが安定して誤差が小さくなります。適切な座標変換は、精度そのものを改善します。

左の棒グラフでは、厳密値・(a)・(b) の3本が目視ではまったく同じ高さで並び、(c) だけが $5.548$ と突出しています。ヤコビアン $r$ を落とすことが「桁を間違える」レベルの誤りであることが一目瞭然です。右は乱数の種を変えて60回試行し、相対誤差の二乗平均平方根を標本数 $M$ に対してプロットしたものです。両者とも $M^{-1/2}$ の破線と平行 — モンテカルロ法の定石どおりの収束率 — ですが、極座標（緑）のほうが常に下にあり、平均で約 $2.5$ 分の1の誤差に収まっています。

実験3：可逆変換で確率密度がどう変わるか

最後に、密度の変換公式(8)を数値検証します。まず1次元の対数正規分布から。$X\sim\mathcal{N}(0,1)$、$Y=e^X$ で、理論密度は $p_Y(y)=\frac{1}{y\sqrt{2\pi}}e^{-(\log y)^2/2}$ でした。サンプルのヒストグラムと重ねます。

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(1)
X = rng.standard_normal(500_000)
Y = np.exp(X)                       # 可逆変換 y = e^x

edges = np.linspace(0.05, 6, 80)
c = 0.5*(edges[:-1] + edges[1:])
# 変換公式：p_Y(y) = p_X(log y) * |dx/dy| = p_X(log y) * (1/y)
p_theory = np.exp(-np.log(c)**2/2) / np.sqrt(2*np.pi) / c

plt.figure(figsize=(10, 4))
plt.subplot(1, 2, 1)
plt.hist(X, bins=80, density=True, color="tab:blue", alpha=0.6, label="標本のヒストグラム")
xs = np.linspace(-4, 4, 400)
plt.plot(xs, np.exp(-xs**2/2)/np.sqrt(2*np.pi), "k-", lw=2, label="理論密度 $p_X$")
plt.title("変換前：標準正規分布 $X$"); plt.xlabel("$x$"); plt.legend(); plt.grid(alpha=0.3)

plt.subplot(1, 2, 2)
plt.hist(Y, bins=edges, density=True, color="tab:orange", alpha=0.6, label="標本のヒストグラム")
plt.plot(c, p_theory, "k-", lw=2, label="変換公式による $p_Y$")
plt.title("変換後：$Y=e^X$（対数正規分布）"); plt.xlabel("$y$"); plt.legend(); plt.grid(alpha=0.3)
plt.tight_layout(); plt.show()

print(f"最大絶対差 = {np.abs(np.histogram(Y, bins=edges, density=True)[0] - p_theory).max():.4f}")
print(f"平均絶対差 = {np.abs(np.histogram(Y, bins=edges, density=True)[0] - p_theory).mean():.5f}")
print(f"理論密度のピーク値 = {p_theory.max():.4f}")

左の図は、下側の $p_X$（左右対称）が写像 $y=e^x$ を通って左側の $p_Y$（右に長い裾）に変わる様子です。$x=-0.8$ 付近の幅 $0.3$ の区間は像では $0.52$ 倍に縮むので密度は $1.91$ 倍に濃くなり、$x=1.1$ 付近では逆に $3.50$ 倍に伸びるので密度は $0.29$ 倍に薄まります。この倍率がちょうど $1/y$ です。右は50万点の標本ヒストグラムと変換公式の重ね描きで、両者は最大絶対差 $0.0299$、平均絶対差 $0.00626$ で一致しています。

出力は「最大絶対差 = 0.0299 / 平均絶対差 = 0.00626 / 理論密度のピーク値 = 0.6567」でした。ピーク値 $0.66$ に対して平均絶対差が $0.006$ ですから、1%程度のずれ — ヒストグラムの標本ゆらぎで説明できる範囲です。左右対称だった正規分布が、右に長い裾を持つ非対称な分布に変わっている 点にも注目してください。指数関数は $x$ が大きいところで急激に引き伸ばす（$dy/dx = e^x$ が大きい）ので、その領域の密度は薄く広がり、逆に $x$ が小さいところは圧縮されて密度が濃くなります。$1/y$ という因子が、まさにこの伸縮を打ち消す形で入っています。

続いて2次元のアフィンカップリング（正規化フローの基本部品そのもの）で確かめます。$\bm{Z}\sim\mathcal{N}(\bm{0},\bm{I}_2)$ から

$$ y_1 = z_1,\qquad y_2 = z_2\,e^{s(z_1)} + t(z_1),\qquad s(u)=0.6\sin u,\ \ t(u)=1.5\sin u $$

と変換します。ヤコビ行列は $\begin{pmatrix}1 & 0\\ \ast & e^{s(z_1)}\end{pmatrix}$ という下三角行列なので $\det\bm{J} = e^{s(z_1)}$、逆写像は $z_1=y_1,\ z_2=(y_2-t(y_1))e^{-s(y_1)}$ です。

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(1)
s = lambda u: 0.6*np.sin(u)      # スケール（対数）
t = lambda u: 1.5*np.sin(u)      # シフト

Z = rng.standard_normal((1_000_000, 2))
Y1 = Z[:, 0]
Y2 = Z[:, 1]*np.exp(s(Z[:, 0])) + t(Z[:, 0])     # アフィンカップリング

# 理論密度：p_Y(y) = p_Z(g^{-1}(y)) * |det J_{g^{-1}}| = p_Z(z) * exp(-s(y1))
xe = np.linspace(-3.5, 3.5, 61); ye = np.linspace(-4.5, 4.5, 81)
xc = 0.5*(xe[:-1] + xe[1:]);     yc = 0.5*(ye[:-1] + ye[1:])
GX, GY = np.meshgrid(xc, yc, indexing="ij")
Z2 = (GY - t(GX))*np.exp(-s(GX))                 # 逆写像
p_theory = np.exp(-(GX**2 + Z2**2)/2)/(2*np.pi) * np.exp(-s(GX))

H, _, _ = np.histogram2d(Y1, Y2, bins=[xe, ye], density=True)

fig, axes = plt.subplots(1, 3, figsize=(14, 4.2))
axes[0].hist2d(Z[:, 0], Z[:, 1], bins=[xe, ye], density=True, cmap="viridis")
axes[0].set_title("変換前：標準2次元正規分布")
axes[1].imshow(H.T, origin="lower", extent=[xe[0], xe[-1], ye[0], ye[-1]],
               aspect="auto", cmap="viridis")
axes[1].set_title("変換後：標本のヒストグラム")
axes[2].imshow(p_theory.T, origin="lower", extent=[xe[0], xe[-1], ye[0], ye[-1]],
               aspect="auto", cmap="viridis")
axes[2].set_title("変換後：変数変換公式による理論密度")
for ax in axes:
    ax.set_xlabel("$y_1$"); ax.set_ylabel("$y_2$")
plt.tight_layout(); plt.show()

print(f"最大絶対差 = {np.abs(H - p_theory).max():.4f}")
print(f"平均絶対差 = {np.abs(H - p_theory).mean():.6f}")
print(f"理論密度のピーク値 = {p_theory.max():.4f}")

左の同心円状の分布が、中央・右のバナナ状に曲がった分布へと変わっています。中央（実際にサンプリングした100万点のヒストグラム）と右（変数変換公式が予言した密度）は、カラースケールを揃えて並べても差が見分けられません。$y_1$ の値に応じて分布の中心が上下にずれるのがシフト $t(y_1)$ の効果、太さが変わるのがスケール $e^{s(y_1)}$ の効果です。

出力は「最大絶対差 = 0.0157 / 平均絶対差 = 0.000551 / 理論密度のピーク値 = 0.1872」となりました。真ん中（実際にサンプリングした分布）と右（変数変換公式が予言した密度）がほぼ同一の絵になっていることが最大の成果です。円形だった等高線が、$y_1$ に応じてバナナのように曲がり、太さも変化しています。 この「曲がり」がシフト $t(y_1)$、「太さの変化」がスケール $e^{s(y_1)}$ の効果で、後者がそのままヤコビ行列式 $\det\bm{J}=e^{s}$ になっています。最大絶対差がピーク値の8%程度あるのは、密度が高い領域でヒストグラムの相対ゆらぎが目立つためで、平均絶対差はピーク値の0.3%に収まっています。

ここで行った計算は、正規化フローが1層分の対数尤度を評価するときにやっていることとまったく同じです。式(11)の $\log|\det\bm{J}_{f}|$ が、この例では $-s(y_1)$ という単なるスカラー関数になっている — 三角ヤコビアンの威力が実感できます。

よくあるつまずき

Q. なぜ絶対値が付くのですか。1変数の置換積分には付いていませんでしたが。 $\det\bm{J}$ が負になるのは向きが反転しているときです。面積・体積は非負の量なので絶対値で符号を捨てます。1変数の置換積分に絶対値が見えないのは、$g$ が減少関数のとき積分の上端と下端が入れ替わり、その「区間の向き」が符号を吸収しているからです。実際、$\int_{g(a)}^{g(b)}$ を常に「小さいほうから大きいほうへ」と書き直せば、$|g’|$ が現れます。

Q. $g$ と $g^{-1}$、どちらのヤコビアンを使うのか毎回迷います。 覚え方はひとつ、「積分変数として新しく使う変数を出発点とする写像」 です。$\int f(\bm{x})d\bm{x}$ を $\bm{u}$ で計算し直すなら、$\bm{u}\mapsto\bm{x}$ の写像のヤコビアンを使います（式(3)）。確率密度で $\bm{y}$ の密度を求めるときは、$\bm{y}$ を入力とする写像 $g^{-1}$ のヤコビアンです（式(8)）。迷ったら、極端な例で検算するのが確実です。$\bm{x}=2\bm{u}$（2倍に拡大）なら面積は4倍になるはずで、$|\det\bm{J}|=4$。$\int f\,d\bm{x} = \int f(2\bm{u})\cdot 4\,d\bm{u}$ が正しい向きです。

Q. 全単射でない変換は使えないのですか。 使えますが、公式に修正が要ります。たとえば $x=u^2$ は $u=\pm\sqrt{x}$ の2価なので、定義域を $u>0$ と $u<0$ に分割して別々に適用し、足し合わせます。一般には「$\bm{x}$ の逆像の個数」で重み付けした面積公式（area formula）に拡張されます。実務上は「区分的に全単射になるよう分割する」で足ります。

Q. $\det\bm{J}=0$ になる点があると公式は破綻しますか。 その点の集合が測度ゼロなら問題ありません。極座標の原点、球座標の $\theta=0,\pi$（極軸上）がその例です。逆に、$\det\bm{J}$ が正の測度を持つ領域でゼロなら、そこでは変換が次元を潰しており、変数変換としては成立しません。

Q. 非正方（$m\neq n$）のヤコビ行列ではどうなりますか。 $\mathbb{R}^m$ から $\mathbb{R}^n$（$mグラム行列式 $\sqrt{\det(\bm{J}^\top\bm{J})}$ を使います。これは曲面積分の面積要素で、$m=n$ のとき $\sqrt{\det(\bm{J}^\top\bm{J})}=|\det\bm{J}|$ となり本記事の公式に一致します。

まとめ

本記事では、重積分の変数変換公式に $|\det\bm{J}|$ が現れる理由を、幾何から確率まで一貫して追いかけました。

• 行列式は体積である。 $2\times2$ の $|ad-bc|$ は平行四辺形の面積、$3\times3$ の行列式はスカラー三重積＝平行六面体の体積。$n$ 次元でも「多重線形・交代・正規化」という体積の公理から行列式が一意に定まる。
• 微分可能写像は局所的に線形。 テイラー展開 $g(\bm{u}_0+\Delta\bm{u})\approx g(\bm{u}_0)+\bm{J}_g\Delta\bm{u}$ により、微小立方体は体積 $|\det\bm{J}_g|h^n$ の平行体に写る。これをリーマン和で足し上げると変数変換公式 $\int_X f\,d\bm{x} = \int_U f(g(\bm{u}))|\det\bm{J}_g|\,d\bm{u}$ が出る。
• 極座標の $r$、球座標の $r^2\sin\theta$ は、いずれも三角関数の基本関係 $\cos^2+\sin^2=1$ を使った数行の計算で導ける。幾何的には「半径が大きいほど弧が長い」「極に近いほど緯線が短い」という直感そのもの。
• ガウス積分 は、2乗して2次元化 → 極座標 → ヤコビアン $r$ が置換を可能にする、という流れで $\sqrt{\pi}$ が出る。余分な因子が問題を解いてくれる好例。
• 確率密度は「質量÷体積」 なので、変換で体積が $|\det\bm{J}_g|$ 倍になれば密度は $1/|\det\bm{J}_g|$ 倍になる。多変量正規分布の $\sqrt{\det\bm{\Sigma}}$ も、2次元正規分布から出るレイリー分布の $r$ も、すべてこの帰結。
• 正規化フロー は $\log p_\theta(\bm{x}) = \log p_Z(f_\theta(\bm{x})) + \log|\det\bm{J}_{f_\theta}(\bm{x})|$ を損失にする。ヤコビ行列を三角行列にする設計（アフィンカップリング）で、$O(n^3)$ の行列式計算が $O(n)$ の総和に落ちる。
• 数値検証では、セル面積比と $|\det\bm{J}|$ の相対誤差が $O(h)$ で消えること、極座標MCが直交座標MCと同じ積分値を（しかも60回試行の二乗平均で約2.5分の1の誤差で）与えること、$r$ を落とすと $3.084$ が $5.548$ という別物になることを確認した。

次のステップとして、以下の記事も参考にしてください。

• ヤコビアンと変数変換（重積分の座標変換公式） — 円筒座標を含む各座標系の計算例
• ガウス積分の導出をわかりやすく解説 — 一般形と正規分布との関係
• 正規化フロー（Normalizing Flow）とは — 本記事の公式を土台にした深層生成モデル
• ディープラーニングにヤコビアンが現れる理由 — 逆伝播・確率変換・感度解析における3つの顔
• 行列式とは？定義と幾何学的な意味 — 行列式そのものをもう一度固める