R: 非矩形領域上での数値積分

投稿者: | 2017年7月28日

数値積分

数値積分は、計算機を使って積分の値を数値的に求める方法です。

例えば、 \[
\int_a^b \int_c^d f(x,y) dxdy
\]
で、\(a,b,c,d,f\) の具体的な値は関数が与えられた上で、その値を計算します。

この例は、\((a,b), (c,d)\) を対角線上の頂点とする長方形上での積分ですが、
一般には、四角い領域以外の積分を計算することもあります。
今回は、その四角い領域以外の積分に焦点を当てていきます。

R の標準的な数値積分

Rには数値積分を行う関数 integrate が標準で入っており、
他にも、拡張的な数値積分を行うパッケージも多数存在します。
それらを使えば、簡単に数値積分を実行することはできますが、
基本は矩形領域、つまり四角い領域での積分を扱っています。

一般の数値積分

ところが、実際に積分する場合は、円上で積分するなど、
矩形領域以外で積分したいことがあります。

つまり、一般には下のような積分を実行します。 \[
\int_R f(x) dx \\
R : g(x) < 0
\]
\(f\) は被積分関数で、\(g\) は積分領域 \(A\) を決定する関数です。
例えば、半径 \(r\) の円上での積分を行う場合は、
\(g(x,y) = x^2 + y^2 – r^2\) となります。

R で用意されている積分がどれも矩形領域に限られていたので、
上のような積分ができる関数を作ってみることにします。

モンテカルロ法

今回は、数値計算法としてモンテカルロ法を用いることにします。
その理由は、積分領域が少々複雑になっても、シンプルな方法で積分が可能だからです。

モンテカルロ法の考え方

◯ 関数の平均値と積分

下の図のような関数 \(f(x)\) の積分を考えてみましょう。

\([a,b]\) 上での積分の値は、図で青く塗りつぶした領域の面積になります。

図の右側には、同じ面積を持った長方形を示しています。
それでは、この長方形の高さ \(\bar{f}\) は何を表しているでしょうか?

この高さは \([a,b]\) 上での \(f(x)\) の平均値を表しています。
これは、砂山の凸凹を平らに均すと、もとの凸凹の平均的な高さになるのと同じです。

長方形の面積は「底の長さ\(\times\)高さ」ですぐに求まるので、
積分の値は \[
\int_a^b f(x) dx = \bar{f} \times (b-a)
\]
で求められることが分かります。

つまり、関数の平均の値が分かれば積分の計算ができるということです。

◯ サンプリングによる関数の平均値の計算

積分の計算は関数の平均値を計算することでできることが分かりました。

モンテカルロ法では、(一様)ランダムに \(x\) をサンプリングして、
\(f(x)\) の平均値を求めます。

\(N\) 個のサンプル点を \(x_1, \cdots x_N\) とすると、
対応した関数の値 \(f(x_1), \cdots, f(x_N)\) が得られます。

そして、その関数の平均値で \(\bar{f}\) を推定するのです。 \[
\bar{f} = \frac{1}{N}\sum_i f(x_i)
\]

サンプル数 \(N\) が十分に大きければ、この値は正しい平均値に近づくので、
この \(\bar{f}\) を使って正しい積分の値を計算することができます。

◯ 非矩形領域での積分

このモンテカルロ法の考え方は多次元、
そして、四角い領域でなくとも同じように使えます。

例えば、下の図のような 2 次元の円上での
関数 \(f(x,y)\) の積分を考えてみましょう。

先程の考え方と同様に、凸凹を均して平らな面をつくると、
右のような円柱に変わります。
このとき、積分の値、つまり体積は同じです。

そして、右の円柱の高さ \(\bar{f}\) は円上での \(f(x,y)\) の平均値ですから、
先程と同様にサンプリングで平均値を計算できます。
後は計算した \(\bar{f}\) の値と、底面積 \(S\) をかけることで堆積が求まります。

◯ 底面積の計算

ところで、非矩形領域の面積 \(S\) はどうやって求めればよいのでしょうか?
円の場合なら手計算できますが、より複雑な領域の場合は困りそうです。

また、領域上で一様ランダムな点をどのように生成すればよいのでしょうか?
矩形領域の場合は、一様分布でそのまま生成できましたが、
円やその他の領域上で偏りなく点を打つにはどうすればよいのでしょうか。

これらの二つの疑問は、乱数の取捨選択をするという簡単な方法で解決できます。

例えば、下の図の青い積分領域を考えます。

青い領域を覆うように矩形領域(外側の四角)を用意します。
そして、四角の中で一様乱数を生成していきます(グレーの丸)。

生成された各乱数 \(x_i\) が青い領域上にあるかどうかは、
\(g(x_i)<0\) であるかどうかを確認すれば分かります。
従って、\(g(x_i)<0\) となった乱数だけを取り出せば、
青い領域上で一様な乱数が生成できます。

この時、青い領域上に載った乱数の割合は、
矩形領域と青い領域の面積の比と等しいはずです。
従って、 \[
\frac{S_青}{S_□} = \frac{N_{青}}{N_□}
\]
から、青い領域の面積を求めることができます。
\(S_青, S_□\)はそれぞれの領域の面積、\(N_青,N_□\)はそれぞれの領域上にある点の数を表します。)

モンテカルロ法の実装

以上に示した方法で数値積分を行う関数のコードは下の通りです。

◯ 関数の引数

  • f : 被積分関数
  • g : 領域を定義する関数
  • ll, ul: g で定義される領域を囲う矩形領域の最小値・最大値
  • N : 生成するサンプル数
  • d : f(x) の変数 x の次元

特に、矩形領域は、d 次元の場合に、
rep(ll,d), rep(ul,d) を対角線上の頂点とする領域となります。

◯ 関数の出力

  • ans : 積分値
  • sd : 積分値の標準偏差

計算実行例

それでは、上の関数を使って数値積分を行ってみます。

例として、次のような積分を扱います。 \[
\int_A xy^2 dxdy \\
A : 1-x-y<0, y\in[0,1], x\in[0,1]
\]
この積分の答えは手で計算することができて、
正確な値は \(\frac{3}{20} = 0.15\) になります。

それでは、サンプル数 10000 で数値積分してみます。
(ただし、サンプルのうち約半数は \(A\) の領域外なので、平均値の計算には使われません。)

積分の結果 ans は正解である 0.15 に近い値になっています。
より精度を高めるには、サンプル数を増やす必要がありますが、 それと同時に計算時間も増えることになります。

MISER 法

サンプル数を増やす以外に精度を高める方法はあるでしょうか?
実はサンプルを増やさなくても、
サンプルの取り方を工夫することで精度を高めることができます。

ここでは、モンテカルロ法の拡張法の一つである MISER 法を紹介します。

MISER法の考え方

◯ どこでも同じ数のサンプルが必要か?

次のような関数の積分を考えてみます。

左半分は、起伏があり、右半分は平坦になっていることがひと目で分かります。

モンテカルロ法による積分は、凸凹を均すことと説明しましたが、
今の場合、右半分はすでに均されているので、
あまり努力しなくてよさそうですよね。

普通にモンテカルロ法を実行すると、右も左も関係なく、
一様にランダムな点を取って全体の平均を求めます。
右だけの平均なら、サンプル 1 つでも十分求まるのに、
たくさんのサンプルを使って平均を求めるのはサンプルの無駄遣いですよね。

つまり、関数の形によって、
サンプルを必要としている領域と、
サンプルをあまり必要としない領域に差がある
のです。

◯ サンプルの数を分配する

MISER 法では、積分領域を分割し、
サンプルを必要としている側のサンプル数を増やし、
サンプルを必要としない側にサンプル数を減らすことで、
同じサンプル数でも効率よく積分の値を計算できるように、
サンプルの分配を行います。

先程の例では、下の図のように、
右側と左側でサンプル数が調整されます。

つまり、平均が求まりやすい右側はサンプル数が少なく、
平均が求まりやすい左側はサンプル数が多くなるように調整します。

具体的なサンプル数は、積分の値の標準偏差に比例して決められます。
つまり。領域を領域1,領域2 の二つに分割した時、
領域1,2 の積分の標準偏差を\(\sigma_1,\sigma_2\) とすると、
領域1,2 のサンプル数 \(N_1, N_2\)\[
N_1 : N_2 = \sigma_1 : \sigma_2
\]
となるように決定されます。

◯ どの方向に分割するか

変数が1次元の場合は、左右の分割を繰り返すだけですが、
多次元になると、横に分割するのか、縦に分割するのか、奥に分割するのか、
という風に、分割する方向も選択肢として現れます。

分割する方向によっては、あまり意味のない分割になる可能性があります。
例えば、右と左では同じくらい凸凹しているのに左右に分割しても、
結局どちらも同じサンプル数を必要とするので意味がありません。

逆に、奥と手前で凸凹の度合いが違うのであれば、
奥の方向に分割することで、サンプルを再分配できて、
精度が向上します。

このように、分割の方向によって、精度の向上しやすさが異なる可能性があります。

そこで、MISER 法ではあらゆる方向で分割してみて、
精度がどのくらい向上するかを調べます。
そして、最も精度が向上する方向に分割するという方法を取ります。

計算実行例

それでは、MISER 法による改良で
本当に、数値積分の精度が向上するかを計算して調べてみましょう。
(MISER 法のコードは長くなったので、記事の一番下に載せます。)

積分は先程モンテカルロ法の計算例で用いたものと同じもので、 サンプル数も同じく 10000 としています。
分割は再帰的に3回繰り返しています。

積分の正しい値は0.15でしたが、
モンテカルロ法よりもさらに近い値を求められていることが分かります。

実は矩形領域の積分でも非矩形領域の積分ができる

さて、今回は \(g(x) < 0\) で表される
一般の領域上での数値積分を行う方法を見てきました。

ところで、R のパッケージでよく実装されている矩形領域の積分で
同じような積分を計算することはできないのでしょうか?

答えは、できます。

被積分関数に領域判定の役割を持たせる

その方法は、被積分関数 \(f\) を変形して、
\(g\) も組み込んだ新しい被積分関数を導入する方法です。 具体的には、 \[
F(x) = \begin{cases}
f(x) & (g(x) < 0) \\
0 & {\rm otherwise}
\end{cases}
\]
という新しい被積分関数 \(F\) を用いて数値積分を行います。

この \(F\) は、積分領域にある時だけ \(f(x)\) の値を返し、
積分領域の外では 0 の値を返します。

つまり、モンテカルロ法でサンプル点を取捨選択した時と同じ役割を、
被積分関数に持たせた上で矩形領域上で積分するのです。
そうすると、上で説明したものと同様の計算を行うことになるので、
非矩形領域での積分が計算できます。

計算実行例

上のような方針で実際に計算してみましょう。

積分は上で使ってきたものと同じものとします。

数値積分には MISER 法を用いており、 矩形積分は、MISER法の引数を \(g(x)=-1\) とすることで実行しています。

非矩形で計算した場合と同様の結果になっていることが確認できます。

MISER 法のコード

montecarlo_integrate とは異なり、
ll,ul はベクトル(座標)で与えます。

また、引数 depth は分割の回数を表します。
ただし、1度分割すると、それぞれの領域について再帰的に分割を開始するので、
より正確には、再帰の深さを表します。


 

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です