JavaScriptでモンテカルロ積分をしてみる

JavaScriptの勉強がてら、webブラウザ上でモンテカルロ積分を実行するプログラムを書いてみました。

モンテカルロ積分の復習
使ってみよう
JavaScriptによる実装

モンテカルロ積分の復習

区間 $[a, b]$ 上で定義された任意関数 $f(x)$、および区間上の一様な確率密度関数 $p(x) = 1/(b – a)$ を考えます。この時、$f$ の積分
$$I := \int_a^b dx\,f(x) = (b-a)\int_a^b dx\,f(x)p(x)$$
を確率論的な観点から近似的に計算してみましょう。

上記の式から明らかなように、積分 $I$ は $f(x)$ の確率密度 $p(x)$ による期待値に比例します。そこでこの期待値を、確率密度 $p(x)$ に従ってサンプリングされた $M$ 個の標本点 $\{x_1, x_2, \cdots, x_M\}$ で近似するというのが、モンテカルロ積分の基本的なアイディアです。

今 $p(x)$ は一様分布なので、標本点は区間 $[a, b]$ から一様にサンプリングすることで得られます。
$I$ の内、$f$ の期待値に相当する部分を $f(x)$ の標本点による平均で置き換えると
$$I_M:= \frac{(b-a)}{M}\sum_{i=1}^M f(x_i) $$
となります。この時、大数の法則より
$$ \lim_{M\to\infty}I_M= I$$
が保証されるため、$M$ が十分大きければ $I_M$ は $I$ の良い近似であることが期待されるのでした。

使ってみよう

被積分関数として、5次までの多項式を用意しました。
計算ボタンを押すたびに乱数が変わるので、違う答えが出ます。

積分範囲
$a =$ $b =$
標本点の数
$M =$ (min : 1, max : 100000)
被積分関数
$f(x)=$ $x^5 +$ $x^4 +$ $x^3 +$ $x^2 +$ $x +$

JavaScriptによる実装

こんな感じのコードで作れます。

<style>
  input[type="number"] {
    width: 70px;
  }
</style>

//入力部分
<p>
  積分範囲<br>
  $a =$
  <input id="range_a" type="number" value="0" step="any" required>
  $b =$
  <input id="range_b" type="number" value="1" step="any" required>

  <br>標本点の数<br>
  $M =$
  <input id="number_M" type="number" value="100" step="1" min="1" max="100000" style="width: 150px;" required>
     (min : 1, max : 100000)

  <br>被積分関数<br>
  $f(x)=$
  <input id="coeff5" type="number" value="0" step="any" required>
  $x^5 +$
  <input id="coeff4" type="number" value="0" step="any" required>
  $x^4 +$
  <input id="coeff3" type="number" value="0" step="any" required>
  $x^3 +$
  <input id="coeff2" type="number" value="0" step="any" required>
  $x^2 +$
  <input id="coeff1" type="number" value="1" step="any" required>
  $x +$
  <input id="coeff0" type="number" value="0" step="any" required>
  <br>
</p>

//計算実行のボタン
<p><input type="button" value="計算" id="checkButton"></p>

//結果の出力場所
<p id="msg1"></p>
<p id="msg2"></p>

//JavaScript
<script>
  //被積分関数 
  function integrand_f(c, x){
      var fx = 0.0;
      var xk = 1.0;

      for (k = 0; k < 6; k++){
          fx += c[k] * xk;
          xk *= x;
      }
      return fx;
  }

  //モンテカルロ積分
  function MonteCarlo(a, b, c, M){
      var I = 0.0;
      var x;

      for (var k = 0; k < M; k++){
          x = Math.random() * (b - a) + a;
          I += integrand_f(c, x);
      }
      I *= (b - a) / M
      return I;
  }
  
  //厳密解の計算
  function exact_solution(a, b, c){
      var I = 0.0;
      for (var k = 0; k < 6; k++){
          I += c[k] * (b**(k+1) - a**(k+1)) / (k + 1.0);
      }
      return I;
  }
  
  //ボタンを押したときに実行される関数
  function butotnClick(){
      var a = parseFloat(range_a.value);
      var b = parseFloat(range_b.value);
      var M = parseInt(number_M.value, 10);

      var c = [];
      for (var i = 0; i < 6; i++){
          var c_dummy = parseFloat(coeff[i].value);
          c.push(c_dummy);
      }

      var I_MC;
      I_MC = MonteCarlo(a, b, c, M);
      var I_EX;
      I_EX = exact_solution(a, b, c);

      msg1.innerText = 'モンテカルロ積分の結果：' + I_MC;

      msg2.innerText = '厳密解：' + I_EX;
  }
  
  //データの受け渡し
  let range_a = document.getElementById('range_a');
  let range_b = document.getElementById('range_b');
  let number_M = document.getElementById('number_M');

  let coeff =[];
  coeff.push(document.getElementById('coeff0'));
  coeff.push(document.getElementById('coeff1'));
  coeff.push(document.getElementById('coeff2'));
  coeff.push(document.getElementById('coeff3'));
  coeff.push(document.getElementById('coeff4'));
  coeff.push(document.getElementById('coeff5'));
  
  let msg1 = document.getElementById('msg1');
  let msg2 = document.getElementById('msg2');

  let checkButton = document.getElementById('checkButton');
  checkButton.addEventListener('click', butotnClick);
</script>

冗長な部分が色々とあるかもしれませんが、その辺はご容赦ください。