p5.jsで円内の領域の上の均等分布を実装する

はじめに

　均等分布。たとえば三角形の上に点を均等にばらまくとしよう。次のようなコードが考えられる。

function setup() {
  createCanvas(400, 400);
  noStroke();
  fill(255);
  background(0);

  for(let i=0; i<20000; i++){
    const x = random(400);
    const y = random(x);
    circle(x,y,1);
  }
}

0～400の範囲の乱数を取り、取った乱数に対して$0$～$x$の範囲の乱数を取る。自然ですね。結果は：

左にだいぶ偏ってしまった。まあ当然か。だって縦方向の乱数の取り方は均等だけれど、その長さが違う。違うにもかかわらず採取ポイントは均等に選ばれる。結果、長さの違うそれぞれの棒に対してその上の点が均等に選ばれるので、狭いところほど密集してしまうのだ。

　これを何とかするのが、この記事のねらいである。

離散化する

　たとえば次のように考えてみる。

x=1,2,3,...,10

に対して、$1,2,3,...,10$の値を取る関数$f(x)$について、$(x,f(x))$の集合を考えるとき、これをランダムで一つ取るという風に考える。このとき、均等にすべての点が選ばれるようにするにはどうすればいいのか。単純にこの点の数は

1+2+3+\cdots +10=55

なので、$1/55$の確率ですべての点は選ばれるのが理想的である。しかしこれを、ランダムで$1,2,3,..,10$の中から一つ選び、その値$x$に対して$1,2,3,...,x$からランダムに選ぶ、としてしまうと、数が少ないほどバリエーションが少なくなってしまい、上の図のような感じになることがなんとなくわかると思う。特に$(1,1)$はこの場合$1/10$の確率で選ばれるので、全く均等ではない。

　そこでまず、この$55$個の数、という意味で、ランダムに$1,2,3,...,55$から一つ選ぶ、と考えてみる。そしてその数$A$に対して、面積の関数：

S(x)=1+2+3+\cdots +x~~~~~(S(0)=0)

を考えて、$S(x-1)<A\leq S(x)$を満たす$x$を考えてみる。この場合、$x$が$1$になる確率は$1/55$、$2$になる確率は$2/55$、...、$10$になる確率は$10/55$である。そしてこの$x$に対して$1,2,3,...,f(x)$からランダムで点を取り、$(x,f(x))$を作ることを考えると、このときの$(1,1)$が選ばれる確率は正しく$1/55$となっている。また$x$が$10$の場合も、$10$になる確率が$10/55$でそれに対し各々の点の選ばれる確率が$1/10$だから、

\frac{10}{55}\times \frac{1}{10}=\frac{1}{55}

で、やはり$1/55$である。このように、すべての点が均等に選ばれる。

　以上のことから、方針としてはまず、$0$から$x$までの面積を与える関数$S(x)$を作る。ただし$x$は$0$から$R$まで動くとしたとき、$S(R)$が全体の面積となるように、適切に選ぶ。そのうえで$0\leq A\leq S(R)$の範囲でランダムに数値を取る。そして$A$に対し

S(x-dx)<A\leq S(x)

となる$x$...というかこれはもはや$y=S(x)$の逆関数なので、要は逆関数を出してそこに$A$をぶち込めばいい。こういうこと。

y=S(x)~~\Leftrightarrow~~x=T(y),~~~~~~~~x=T(A)

で、その$x$について$0\leq y \leq f(x)$を満たす$y$をランダムに選べばいいと思う。すなわち、

  x=T(random(S(R))), y=random(f(x))

とすると、$(x,f(x))$は均等に選ばれる。$S(R)$は要するに総面積である。$T$は面積関数の逆関数だが、このような関数は基本、積分で定義される。

S(x)=\int_0^x f(t)dt.~~~~~y=S(x),~~x=T(y).

だから積分で面積関数を作って逆関数を取ればいいですね。

実装

　というわけで三角形の問題に戻る。$x$は$0$から$400$まで動く。面積関数は、

S(x)=\int_0^x tdt = \frac{x^2}{2},~~~x=\sqrt{2y}

なので$\sqrt{2y}$である。総面積は$400\times 400\times 0.5=80000$なので、要するに$160000$に通常の小数ランダムを掛けて平方根を取れば望みの$x$が得られることになる。この$x$に対して縦の棒から均等に点を取ればいい。今の場合は単純にこの$x$に小数ランダムを掛ければいい。実装すると、こう：

  for(let i=0; i<20000; i++){
    const x = sqrt(random(160000));
    const y = random(x);
    circle(x,y,1);
  }

実行結果：

偏りが消えましたね。

　円の場合もやってみる。半径200の円の中の点を均等に取ってみよう。まず単純に0～200のどれかを取り、中心からの距離とし、それとは別に角度を0～$2\pi$でランダムに取る。こんな風に：

  for(let i=0; i<20000; i++){
    const r = random(200);
    const t = random(TAU);
    circle(200+r*cos(t),200+r*sin(t),1);
  }

結果は：

見事に中心によりましたね。理由は同じ。そこでさっきの理論を使う。まずこの場合$x$は0から200まで動く。それに対する面積関数は簡単に出る。

S(x)=\pi x^2,~~~~x=T(y)=\sqrt{\frac{y}{\pi}}

逆関数も簡単に出る。そして面積の最大は$40000\pi$なので、愚直にこうすればいい。今回は$x$を出した場合のスライスされた領域はすべて円環なので、$x$に依らずランダムで角度を取得するだけでいい。こんな感じ：

  for(let i=0; i<20000; i++){
    const r = sqrt(random(40000*PI)/PI);
    const t = random(TAU);
    circle(200+r*cos(t),200+r*sin(t),1);
  }

もっとも$\pi$を相殺するなどいろいろ無駄が多いが、説明のためである。あんま省略すると分かりにくくなるので。それで、この場合の実行結果はこちら：

さっきと全然違いますね。そういうことですね。

　要は$x$での断面から均等に点を取ればいい。球面上の場合はこんな感じになる。

  for(let i=0; i<4000; i++){
    const s = acos(1-random(4*PI)/TAU);
    const t = random(TAU);
    positions.push(createVector(sin(s)*cos(t),sin(s)*sin(t),cos(s)).mult(200));
  }

これで半径$200$の球面上の点をランダムに取得できる。実際にこれで3D空間に球を並べたりできるが、p5の通常のwebGLの枠組みでこれをやると激重になってしまうので、インスタンシングをつかうなど、何かしら工夫が必要かもしれないですね。

おわりに

　カヤックさんの次の記事は非常に参考になります。一度読んでみるといいです。面白いです。
　円や球面にランダムに点を散らしたい〜高校数学を使って分布関数を作る〜

　ここまでお読みいただいてありがとうございました。

補足

　上の記事でも扱っているが、面積計算の際に重みを付けることで、分布を変えることができる。その場合、ランダム値の上限は面積ではなく全体の積分値となるのでそこは注意が必要である。たとえば円の分布について、

  for(let i=0; i<20000; i++){
    const r = 200*pow(random(1),1/6); // ここの乗数を1/2ではなく1/3,1/4,...としてみる
    const t = random(TAU);
    circle(200+r*cos(t),200+r*sin(t),1);
  }

などとしてみる。実は1/2だと均等になるが、1/6ではどうなるかというと、

かなり端っこに寄る。逆に2や3にすると中心に寄る。いろいろ遊んでみると楽しいです。