Wilson 分布を確認していたら、「確率変数の和の小数部分の分布」やら「二次元ランダムウォークの成分の独立性」が問題となった

原子が単位格子中にランダムに散らばっているとしよう。i番目の原子の座標(fractional coordinate)を $\mathbf{x}_i$ 、原子散乱因子を $f_i$ とすると、結晶全体からの散乱は指数を h として、

$F(\mathbf{h}) = \sum_i f_i e^{2\pi i \mathbf{h} \cdot \mathbf{x}_i}$

で表される。なお、この手の式ではいちいち書かないのが普通だが、 $f_i$ は指数 h に依存するから、正確には $f_i(\mathbf{h})$ という意味である。

もっとも簡単に、原子がすべて同種である場合のある指数の反射をシミュレートしてみよう。この場合、 $f_i$ は定数なので、1としてしまっても構わない。h と x の内積を取っている部分は、x の成分が[0, 1)の一様分布であるが、その3次元分の和なので、中心極限定理により正規分布に近づく。しかし、 $2\pi i$ をかけて位相項として使うだけなので、和の小数部分だけが問題となって、それはやはり[0, 1)の一様分布となる……ようだ。iid な確率変数の和の小数部分は一様分布に収束するという定理があるらしい(fractional part, sum of random variables などで検索すると、Encyclopedia of Mathematicsなどが見つかる)。今回は、x, y, z 成分ごとに指数で重みがついているので iid とは言えないが、元が一様分布なので、「一様分布の定数倍の小数部分も一様分布」を認めれば、これでよいことになる。＃このあたりの証明は読もうと思ったが、数学的に難しかったので保留

Ntry <- 10000
ret <- rep(0, Ntry)
h <- c(10, 3, 4) # 指数。ここをどんな整数の組にしても...

for (i in 1:Ntry) {
  ret[i] <- sum(runif(3, min=0, max=1) * h) %% 1
}
hist(ret) # 一様分布!

結局、全体の構造因子の分布は、二次元ランダムウォーク(歩幅は1で固定だが、方向は一様分布)をNatom歩したあとの到着点の分布と同じになる。R の複素数型は賢いので、かなりコンパクトに書ける。

Natom <- 1000
Ntry <- 10000

ret <- rep(0, Ntry)
for (i in 1:Ntry) {
   theta <- runif(Natom, min=0, max=2*pi)
   ret[i] <- sum(complex(arg=theta))
}
hist(Re(ret))
plot(ret, cex=0.1)

原点を中心とする二次元正規分布になっているように見える。

まず実部について考えよう。1つ1つの原子からの寄与は、-1 から 1 まで、 $\cos \theta$ にしたがって分布する確率変数である。その平均は、 $\int_0^{2\pi} {1 \over {2\pi}} \cos \theta d\theta = 0$ であり、分散は、 $\int_0^{2\pi} {1 \over {2\pi}} \cos^2 \theta d\theta = {1 \over 2}$ である。この分布は正規分布ではない。それどころか、両端 -1 と 1 で大きく、中央 0 で最小値となる変な分布である(確率変数の変数変換 - biochem_fan's note で実験した)。

結晶全体の構造因子は、この確率変数 Natom個分だから、平均と分散の線形性により、平均が 0、分散が $N_{atom} \over 2$ となる。Natom が十分に大きい場合は、中心極限定理によってこの分布は正規分布に漸近する。まったく、驚くべきことではないか!

虚部についても同様。

この後、「実部と虚部は独立なので」といって共分散成分のない二次元正規分布を持ち出すのが教科書でおきまりのストーリなのだが、本当に独立なのかは自明ではない。少なくとも、原子1つ1つの寄与は $\cos \theta$ 対 $\sin \theta$ という強い関係性で結ばれているからである。一つの直感的な説明は、片方がある値だと分かっても、もう片方に正負の自由度があるため、ちょうど相関が消えてしまうというものである。だが、それでいいのだろうか。＃というところで、未解決