読者です 読者をやめる 読者になる 読者になる

2/2 複素標準正規分布と複素標準ガウス雑音について

複素標準ガウス雑音は、実軸方向と虚軸方向にそれぞれ分散0.5のガウス雑音を足せば作れる。
ここでの標準は、平均0で分散1を意味する。
これを示すために、実軸方向と虚軸方向にそれぞれ分散0.5のガウス雑音を生成し、(正規分布に従う乱数がすなわちガウス雑音であることに留意)その同時確率分布をヒストグラムで近似して、複素標準正規分布の解析的な式と比較する。

分散なので{ \displaystyle
\sigma^2 = 0.5
}
であることに留意。

式はzを複素数として、
{ \displaystyle
p(z) = \frac{1}{\pi}e^{-|z^2|}
}
で与えられる。
サンプルコード。

import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import pylab as pl
from pandas import *
size = 1000000

#乱数生成
x = np.random.normal(0.0,np.sqrt(0.5),size)
y = np.random.normal(0.0,np.sqrt(0.5),size)

#描画用のメッシュを切っている
xmesh = np.arange(-3.5,3.5,0.1)
ymesh = np.arange(-3.5,3.5,0.1)

BIN = len(xmesh)

#ヒストグラム作成。pl.histが二次元配列を返すのはあまり触れられていない。調べればわかるからか?なんにせよ便利だ。
xhist = pl.hist(x,bins = BIN,range=(-3.5,3.5),normed = True)
yhist = pl.hist(y,bins = BIN,range=(-3.5,3.5),normed = True)


Z = np.zeros((BIN,BIN))

#ノイズの分布をヒストグラムで近似している。真のpython使いなら頑張って一個forループを削る気がする。
for i in range(BIN):
    for j in range(BIN):
        Z[i][j] = xhist[0][i] * yhist[0][j]

#これで描画用の配列を作るらしい。よくわからん。
#>やすうみさんから指摘:ZがBIN*BINのサイズをもつ配列なので、それと一対一対応で点をプロットするのにBIN*BINサイズの配列を作成している。
X,Y = np.meshgrid(xmesh,ymesh)

Z2 = np.zeros((BIN,BIN))

#複素標準正規分布を直接計算。やはり真のpython使いはforループを削るはず。わたしはそんなことするならCで書く。
for i in range(BIN):
    for j in range (BIN):
        Z2[i][j] = 1.0/np.pi*np.exp(-xmesh[i]*xmesh[i] - ymesh[j]*ymesh[j])

#グラフ作成
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_wireframe(X,Y,Z,color ='red')
ax.plot_wireframe(X,Y,Z2,color = 'blue')
plt.show()

グラフはこうなる。

f:id:glamorousdammy:20150202033425p:plain

教訓:目の前に計算機環境があるのにこういったことを調べないのは怠慢である、おかげで効率が落ちたことを自覚し、きちんと調べつくすべきだ、今後は気をつけねばならない。