点を打つのはいいのですが、本当にそれが超球面上に一様に分布しているのか確認したいところです。
ある点を固定して、超球の中心を挟んでその点と他の点との成す角度がある分布をしているはずです。例えば3次元で考えると、角度θとθ+dθの間の表面積は
2πr2sinθdθ
（rは半径）となるから、
f(θ) = sinθ/2
という分布関数になるはずです。一般の次元Dで超表面積は、
V_{D-2}(rsinθ)dθ
（V_D(r)は半径rの超球の体積）だから、
sin^{D-2}θ
に比例するはずです。
D = 6なら、
(8/3π)sinθ
です。
結果的には、次のようになりました。理想の分布との比を求めていて、coarseが前回ダメと書いた方法、fineが前回求めた方法です。いちおうcoarseはダメらしいことがわかります。
from itertools import combinations from math import sqrt, log, cos, sin, pi, acos import random def histogram(xs, first, last, M): h = [ 0 ] * M for x in xs: k = int((x - first) / (last - first) * M) h[k] += 1 return h def Box_Muller(): r1 = random.random() r2 = random.random() return sqrt(-2.0 * log(r1)) * cos(2.0 * pi * r2) def random_hyper_sphere(D): xs = [ Box_Muller() for _ in range(D) ] r = sqrt(sum(x * x for x in xs)) return [ x / r for x in xs ] def random_hyper_sphere_coarse(D): xs = [ random.uniform(-1, 1) for _ in range(D) ] r = sqrt(sum(x * x for x in xs)) return [ x / r for x in xs ] def angle(pt1, pt2): return acos(sum(x * y for x, y in zip(pt1, pt2))) # D : 6 def ideal_distribution(M): xs = [ pi * k / M for k in range(M + 1) ] ys = [ 0.375 * x - sin(2*x) / 4 + sin(4*x) / 32 for x in xs ] return [ (y2 - y1) / (0.375 * pi) for y1, y2 in zip(ys, ys[1:]) ] def distribution(points_generator, N, M): pts = [ points_generator(D) for _ in range(N) ] L = N * (N - 1) / 2 angles = [ angle(pt1, pt2) for pt1, pt2 in combinations(pts, 2) ] hist = histogram(angles, 0.0, pi, 10) return [ f / float(L) for f in hist ] D = 6 # dimension N = 1000 # number of samples M = 10 # number of bins dist = distribution(random_hyper_sphere, N, M) dist_coarse = distribution(random_hyper_sphere_coarse, N, M) dist_ideal = ideal_distribution(M) angles = [ pi * (k + 0.5) / M for k in range(M) ] for theta, d, dc, di in zip(angles, dist, dist_coarse, dist_ideal): print "%.6f\t%.6f\t%.6f" % (theta, d / di, dc / di)
