ガウス積分のメモ

ガウス波束や正規分布の取り扱いで何かと出てくる積分です。
さらに多項式が掛かっていると多少やっかいです。
出てくるたびに毎回検索したり調べたり計算したりは面倒なので、一度どういう形になるかをまとめてみようと思います。
ここでは1変数のものを取り扱います。多変数の拡張は比較的簡単かと思います。

今回は数式のみですが、数式が妥当かどうかの確認を裏でjuliaを使って数値計算でやっているので、そのうちそれも記事にしようかなと思います。
追記：しました。https://qiita.com/sage-git/items/07208f54af5c57663ed6

あと、数式の計算はいきなり記事で書きながらLaTeX表記上でやっていて、行間を適切に削るのも面倒なので計算途中も多くを残しています。

状況

中心 $r$、幅 $a^{-1}$ のガウス関数に $x^n$ を掛けて $-\infty$ から $\infty$ まで積分します。 ガウス関数と多項式の中心がずれているので、奇数次のときは奇関数なので０、というわけにはいきません。

基本

$$
I_0(r, a) = \int_{-\infty}^\infty \exp\left(-a(x-r)^2 \right)dx = \sqrt{\frac{\pi}{a}}
$$

よく分子分母がどっちか忘れるのですが、最近 $dx$ と $a$ の次元を考えたら決まることに気づきました。

多項式を掛けて

$$
I_n(r, a) = \int_{-\infty}^\infty x^n\exp\left(-a(x-r)^2 \right)dx
$$

基本の式より$I_0(r, a) = \sqrt{\pi/a}$ です。

すぐに計算できるところ

1次

計算

まず $I_0(r, a)$ を $r$ で微分します。

$$
\frac{\partial }{\partial r}I_0(r, a) = \frac{\partial }{\partial r}\sqrt{\frac{\pi}{a}} = 0
$$

一方で積分と微分を入れかえると

$$
\int_{-\infty}^\infty \frac{\partial }{\partial r}\exp\left(-a(x-r)^2 \right)dx = 2a\int_{-\infty}^\infty (x-r)\exp\left(-a(x-r)^2 \right)dx
$$

右辺のカッコ内を展開すると

$$
2a\int_{-\infty}^\infty x\exp\left(-a(x-r)^2 \right)dx - 2ar\int_{-\infty}^\infty \exp\left(-a(x-r)^2 \right)dx =2aI_1(r, a) - 2arI_0(r, a) = 0
$$

まとめると

$$
I_1(r,a) = rI_0(r,a) = r\sqrt{\frac{\pi}{a}}
$$

応用

$I_1(r, a)/I_0(r, a) = r$ は平均 $r$、分散 $2/a$ の正規分布 $x$ の期待値 に相当します。 $I_0(r, a)$ で割るのは規格化因子です。これが $r$ になるのは直感的にも合っていますね。

2次

計算

今度は $a$ で微分してみます。

$$
\frac{\partial }{\partial a}I_0(r, a) = \frac{\partial }{\partial a}\sqrt{\frac{\pi}{a}} = -\frac{1}{2a}\sqrt{\frac{\pi}{a}}
$$

一方で

$$
\int_{-\infty}^\infty \frac{\partial }{\partial a}\exp\left(-a(x-r)^2 \right)dx = - \int_{-\infty}^\infty (x-r)^2\exp\left(-a(x-r)^2 \right)dx
$$

右辺カッコ内を展開すると

$$
-I_2(r, a) + 2rI_1(r, a) - r^2I_0(r, a) = -\frac{1}{2a}\sqrt{\frac{\pi}{a}}
$$

$$
I_2(r, a) = \frac{1}{2a}\sqrt{\frac{\pi}{a}}+ 2rI_1(r, a) - r^2I_0(r, a) = \left(\frac{1}{2a} + 2r^2 -r^2\right) \sqrt{\frac{\pi}{a}}
$$

まとめると

$$
I_2(r, a) = \left(\frac{1}{2a} + r^2\right) \sqrt{\frac{\pi}{a}}
$$

別解

ちなみに $I_1(r, a)$ を $r$ で微分してみますと

$$
\frac{\partial }{\partial r}I_1(r, a) = \frac{\partial }{\partial r}r\sqrt{\frac{\pi}{a}} = \sqrt{\frac{\pi}{a}}
$$

一方で積分の方は

$$
\int_{-\infty}^\infty \frac{\partial }{\partial r}x\exp\left(-a(x-r)^2 \right)dx = 2a\int_{-\infty}^\infty x(x-r)\exp\left(-a(x-r)^2 \right)dx
$$

$$
2aI_2(r, a) - 2arI_1(r, a)= \sqrt{\frac{\pi}{a}}
$$

まとめると
$$
I_2(r, a) = \left(\frac{1}{2a} + r^2\right) \sqrt{\frac{\pi}{a}}
$$

で、無事 $a$ で微分したときと同じ結果になりました。

応用

$x^2 - r^2$ をガウス関数$\exp\left(-a(x - r)^2\right)$に掛けて積分すると $I_2(r, a) - r^2I_0(r, a) = \frac{1}{2a} I_0(r, a)$ ですが、$x$ を標本平均と思うと標本平均の分散となります（多変数バージョンにしても確認できます）。これは不偏分散を求める際、差の2乗を $N$ ではなく $N-1$ で割ることにつながります。

-1次

計算

普通のガウス関数のフーリエ変換を考えます。$\sqrt{2\pi}$で割らないタイプ。

$$
\int_{-\infty}^\infty \exp\left(-a(x-r)^2\right) \exp(-ikx) dx
$$

指数関数の肩を $x$ に関してまとめなおしますと

$$
-a(x-r)^2 - ikx = -a\left(x^2 - \left(2r - \frac{ik}{a}\right) x + r^2\right)\\ = -a\left(x - \left(r - \frac{ik} {2a}\right)\right)^2 - ikr - \frac{k^2}{4a}
$$

中心位置が虚数なのに目をつぶって積分すると、

$$
\int_{-\infty}^\infty \exp\left(-a(x-r)^2\right) \exp(-ikx) dx = \sqrt{\frac{\pi}{a}}\exp(-ikr)\exp\left(-\frac{k^2}{4a}\right)
$$

普通に中心位置がシフトしたガウス関数のフーリエ変換です。

あと、$x^{-1}$ のフーリエ変換

$$
\int_{-\infty}^\infty \frac{1}{x} \exp(-ikx) dx = - \pi i \mathrm{sgn}(k)
$$

も準備しておきます。証明は省略します。

さらに畳み込みに関するフーリエ変換の性質も確認しておきます。$\hat{f}$ というように $\hat{\ }$ を付けると $f$ のフーリエ変換と表しますと、関数の畳み込みのフーリエ変換は

$$
\int_{-\infty}^\infty \left( \int_{-\infty}^\infty f(y)g(x - y)dy\right)\exp(-ikx)dx = \hat{f}(k)\hat{g}(k)
$$

とフーリエ変換の積で計算できますが、逆にフーリエ変換の畳み込みは関数の積を使って

$$
\int_{-\infty}^\infty f(x)g(x)\exp(-ikx)dx = \frac{1}{2\pi} \int_{-\infty}^\infty \hat{f}(\kappa)\hat{g}(k - \kappa)d\kappa
$$

とできます。 $f(x) = \exp\left(-a(x-r)^2\right)$ と $g(x) = 1/x$ を代入すると

$$
\int_{-\infty}^\infty \frac{1}{x}\exp\left(-a(x-r)^2\right) \exp(-ikx) dx = -\frac{i}{2}\sqrt{\frac{\pi}{a}} \int_{-\infty}^\infty\exp(-i\kappa r)\exp\left(-\frac{\kappa^2}{4a}\right)\mathrm{sgn}(k - \kappa)d\kappa
$$

となります。さらに $k = 0$ とすると、左辺は求める積分になることが分かります。ということで右辺の計算を進めます。積分変数を途中で $\kappa$ から $k$ に変えて

$$
\int_{-\infty}^\infty \frac{1}{x}\exp\left(-a(x-r)^2\right) dx =\frac{i}{2}\sqrt{\frac{\pi}{a}} \int_{-\infty}^\infty\exp(-i\kappa r)\exp\left(-\frac{\kappa^2}{4a}\right)\mathrm{sgn}(\kappa)d\kappa　\\
= \frac{i}{2}\sqrt{\frac{\pi}{a}} \left[-\int_{-\infty}^0\exp(-ikr)\exp\left(-\frac{k^2}{4a}\right)dk + \int_0^\infty\exp(-ikr)\exp\left(-\frac{k^2}{4a}\right)dk\right]\\
=\frac{i}{2}\sqrt{\frac{\pi}{a}} \left[-\int_0^\infty\exp(ikr)\exp\left(-\frac{k^2}{4a}\right)dk + \int_0^\infty\exp(-ikr)\exp\left(-\frac{k^2}{4a}\right)dk\right]\\
=\frac{i}{2}\sqrt{\frac{\pi}{a}} \left[-2i\int_0^\infty\sin(kr)\exp\left(-\frac{k^2}{4a}\right)dk\right]\\
=\sqrt{\frac{\pi}{a}} \left[\int_0^\infty\sin(kr)\exp\left(-\frac{k^2}{4a}\right)dk\right]
$$

と積分できそうな形になりました。係数を移してもうちょっと頑張ります。$\sin(x)$ は指数表記の虚部 $\Im[\exp(ix)]$ ということで、

$$
\sqrt{\frac{a}{\pi}}\int_{-\infty}^\infty \frac{1}{x}\exp\left(-a(x-r)^2\right) dx = \Im\left[\int_0^\infty\exp(ikr)\exp\left(-\frac{k^2}{4a}\right)dk\right] \\
= \Im\left[\int_0^\infty\exp\left(-\frac{1}{4a}\left(k - 2iar\right)^2 - ar^2\right)dk\right] \\
= \exp\left(-ar^2\right)\Im\left[\int_{-2iar}^{\infty - 2iar}\exp\left(-\frac{k^2}{4a}\right)dk\right] \\= 2\sqrt{a}\exp\left(-ar^2\right)\Im\left[\int_{-i\sqrt
{a}r}^{\infty - i\sqrt{a}r}\exp\left(-k^2\right)dk\right]
$$

$\int_{-2iar}^{\infty - 2iar}$ の積分は $x$軸に平行に、$x$軸から $-2iar$ だけ上下にずれたところをたどっていってます。つまり複素平面上での積分になります。このとき被積分関数は性質がよさそうなので、始点終点が同じならば結果は同じだろうということで、実軸をたどるようにします。

$$
\int_{-i\sqrt{a}r}^{\infty - i\sqrt{a}r}\exp\left(-k^2\right)dk = \left(\int_{-i\sqrt{a}r}^{0} + \int_0^{\infty}+\int_{\infty}^{\infty - i\sqrt{a}r}\right)\exp\left(-k^2\right)dk \\ =
\int_0^{i\sqrt{a}r}\exp\left(-k^2\right)dk + \int_0^{\infty}\exp\left(-k^2\right) dk \\=
\frac{\sqrt{\pi}}{2}\mathrm{erf}(i\sqrt{a}r) + \frac{\sqrt{\pi}}{2}
$$

実軸上無限遠でちょっと虚軸を動いた程度の積分は０としました。$\mathrm{erf}(x)$ は普通の誤差関数です。よって

$$
2\sqrt{a}\exp\left(-ar^2\right)\Im\left[\int_{-i\sqrt
{a}r}^{\infty - i\sqrt{a}r}\exp\left(-k^2\right)dk\right]\\
= \exp\left(-ar^2\right)\Im \left[\sqrt{a\pi}\left(1 + \mathrm{erf}(i \sqrt{a} r)) \right) \right] \\
= \exp\left(-ar^2\right)\sqrt{a\pi}\Im \left[\mathrm{erf}(i \sqrt{a} r)\right]
$$

とできます。虚部をとるので $1$ は消えます。元の式に入れてまとめると

$$
I_{-1}(r, a) = \int_{-\infty}^\infty \frac{1}{x}\exp\left(-a(x-r)^2\right) dx = \pi\exp\left(-ar^2\right)\Im \left[\mathrm{erf}(i \sqrt{a} r)\right]
$$

ですが、虚数誤差関数 $\mathrm{erfi}(x) = -i\mathrm{erf}(ix)$を用いると

$$
I_{-1}(r, a) = \pi\exp\left(-ar^2\right)\mathrm{erfi}\left(\sqrt{a} r\right)
$$

となります。正直ぱっと見発散する積分で、本当か？という気分になりますが、これを次回確認してみたいと思います。

追記

$1/x$ のフーリエ変換と畳み込みの話は、投稿当初は間違えていましたが、 @septcolor 様に指摘いただきました。
読みやすさのため修正後のみの内容ですが、間違いのあらすじは $1/x$ のフーリエ変換の係数を $2\pi$ で余計に割っていて、一方関数の積のフーリエ変換で $2\pi$ を忘れていいたため、打ち消しあって結果は正しかったというものです。

応用

この積分は基底関数としてガウス関数系にして、クーロン相互作用を計算したときに似たような形が出てきます。量子化学計算の実装で誤差関数がでてくるのはこれのせいです。
この積分が発散せず値を持つだろうという気持ちは、このことに基づいています。

ただ、クーロン積分は三次元中で、空間中の点から原子核までの距離の逆数 $1/|\vec{R} - \vec{r}|$ のフーリエ変換は $1/k^2$ になって若干違う形（誤差関数が実数で済む）になります。

-2次以降

裏で確かめたのですが、数値積分を $[-\infty, -\epsilon]$、$[\epsilon, \infty]$ に分けて $\epsilon$ を小さくしていくと、-2次と-3次の場合積分の結果がどんどん大きくなりました。(-1次の時は求めた式と同じ値に落ち着きました)

そもそも収束するかどうか怪しいですね。もしかしたら $n < -1$ では $x = 0$ 近傍が大きくなりすぎて発散するのかもしれません。

一般化

というわけで0以上の時だけを考えます。まず漸化式を求めておきます。

$$
I_n(r,a) = \int_{-\infty}^\infty x^n\exp\left(-a(x - r)^2\right)dx
$$

いつも通り微分しますと

$$
\frac{\partial I_n(r, a)}{\partial r} = \int_{-\infty}^\infty 2a\left(x^{n+1} - rx^n\right)\exp\left(-a(x - r)^2\right)dx = 2aI_{n+1}(r, a) - 2arI_n(r, a)
$$

$$
I_{n+1}(r, a) = rI_n(r, a) + \frac{1}{2a}\frac{\partial I_n(r, a)}{\partial r}
$$

$I_0(r, a) = \sqrt{\pi/a}$ なので、どうも $r$ の多項式になりそうです
部分積分でも漸化式を得るのですが、今度は微分がないものの$n+2$, $n+1$, $n$ が絡む漸化式となり、同値かどうかは確かめてないです。

いくつか具体例を挙げると、$I_1(r,a) = rI_0$ と $I_2(r, a) = (r^2 + (2a)^{-1})I_0$ は求めたとおりですが、

$$
I_3 = r\left(r^2 + \frac{1}{2a}\right)I_0 + a^{-1}rI_0 = \left(r^3 +\frac{3}{2a}r\right)I_0 \\
I_4 = r\left(r^3 +\frac{3}{2a}r\right)I_0 + a^{-1}\left(\frac{3}{2}r^2 +\frac{3}{4a}\right)I_0 = \left(r^4 +3a^{-1}r^2+\frac{3}{4}a^{-2}\right)I_0 \\
I_5 = r \left(r^4 +3a^{-1}r^2+\frac{3}{4}a^{-2}\right)I_0 + a^{-1}\left(2r^3 +3a^{-1}r\right)I_0 = \left(r^5 + 5a^{-1}r^3+ \frac{15}{4}a^{-2}r\right)I_0
$$

と続きます。 一般項を思いついた、またはどこかで見つけたら追記したいと思います。

みつけた

$$
\int_{-\infty}^\infty x^n\exp(-px^2 + 2qx)dx = n! \exp(q^2/p)\sqrt{\frac{\pi}{p}}\left(\frac{q}{p}\right)^n\sum_{k=0}^{\lfloor n/2 \rfloor} \frac{1}{(n-2k)!k!}\left(\frac{p}{4q^2}\right)^k
$$

というのを見つけました。¹

$$
-a(x - r)^2 = -ax^2 + 2arx - ar^2
$$

より、 $p = a$、$q = ar$ だといえるので、$ \exp(ar^2)$ をキャンセルして

$$
I_n(r, a) = n!\sqrt{\frac{\pi}{a}}r^n\sum_{k=0}^{\lfloor n/2 \rfloor} \frac{1}{(n-2k)!k!}\left(\frac{1}{4ar^2}\right)^k
$$

です。数値計算チェックも大丈夫でした。

追記

記事冒頭にも書きましたが、数値計算(julia)によるチェックもしました。
https://qiita.com/sage-git/items/07208f54af5c57663ed6

1. Computation of Overlap, Kinetic and Nuclear Attraction Integrals by a Novel Method. Turk J Phy 25 (2001) , 237 – 248 ↩