2024年2月3日土曜日

円の長さ

正方形の長さからの続き

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

練習として,正方形内のランダムな2点の距離の期待値$\ d \ $を計算できることを確認した。次にトライするのが円であるが,ネットで検索しても生成AIにきいてもあまり適切な解答が得られない。一番近いのが,Yahoo知恵袋の「半径1の円内の任意の2点間の距離の期待値は?」だ。これも結局解析的な答えがでなくて,数値計算で $\ 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に投げてみたが,忍耐可能時間内には答えが出なかった。

0 件のコメント:

コメントを投稿