円の長さ

正方形の長さからの続き

ある図形の大きさの指標となる長さを，図形内に一様分布する2点の距離の期待値として定義することで，都道府県の形や大きさを，面積や周長だけでなく"長さ"で特徴づけるという話をしている。

練習として，正方形内のランダムな2点の距離の期待値$\ d \ $を計算できることを確認した。次にトライするのが円であるが，ネットで検索しても生成AIにきいてもあまり適切な解答が得られない。一番近いのが，Yahoo知恵袋の「半径１の円内の任意の２点間の距離の期待値は？」だ。これも結局解析的な答えがでなくて，数値計算で $\ d=0.9054\ $という値を出している。

$d = \int_0^{2\pi} (\frac{1}{2\pi}) d\theta_1 \int_0^{2\pi} ( \frac{1}{2\pi}) d\theta_2 \int_0^1 (2 r_1) dr_1 \int_0^1 (2 r_2) dr_2 \sqrt{r_1^2+r_2^2-2r_1 r_2 \cos(\theta_1-\theta_2) }$

ただし，各積分の()内がそれぞれの変数に対応する確率密度関数，$p(\theta_1)$，$ p(\theta_2)$，$q(r_1)$， $q(r_2)\ $であり，それぞれの変数で積分すると1になるように規格化されている。

これを計算するためには，与えられた4変数の確率分布関数から変数変換によって，積分可能な形に持ち込む必要があるが，なかなか難渋する。しかたがないので，とりあえずJuliaとMathematicaで数値計算してみる。

a=zeros(Float64,1000001,2)

function ju(a,n)

  k = 0

  for i in 1:n

      x = 2*rand()-1

      y = 2*rand()-1

      if x^2+y^2 < 1.0

        k = k + 1

        a[k,1] = x

        a[k,2] = y

      end

  end 

  return k

end

function su(a,n)

  m = div(n,4)*3

  sum = 0

  for i = 1:m

    for j = i:m

      sum = sum + sqrt((a[i,1]-a[j,1])^2+(a[i,2]-a[j,2])^2)

    end

  end

  return sum/binomial(m,2)

end

n=300000

@time su(a,n)

235821
24.223278 seconds (7.18 k allocations: 500.875 KiB, 0.07% compilation time)
0.9055523

生成AIの2つが答えた，$\dfrac{4}{\pi}=1.27324$はたぶん誤っていたということだろう。

解析的に計算できないかと思うのだが，角度積分が完全楕円積分の形になるので，これをさらに積分するのはちょっと難しそうだった。角度積分を後回しにしてもさらに面倒か。Mathematicaに投げてみたが，忍耐可能時間内には答えが出なかった。