1

以下にリストされている出力のコードと写真がありますが、プロットされた特定の標準偏差内でこれらの球からランダムなサンプルを取得したいと考えています。変数 sdwith は、ワイヤーメッシュの出力用のコードでこれを設定するために使用されます。random.multivariate_normal はサンプリングを行いますが、サンプリングする最大確率または標準偏差の数を設定することはできません。これはnumpyで可能ですか、それともこれを行う最善の方法は何ですか?

def sphere(r=1.0,npts=(20,20)):
     """Create a simple sphere.
    Returns x, y, z coordinates for the sphere
    """
     phi=linspace(0,pi,npts[0])
     theta=linspace(0,2*pi,npts[1])
     phi, theta = meshgrid(phi,theta)
     x = r*sin(phi)*cos(theta)
     y = r*sin(phi)*sin(theta)
     z = r*cos(phi)
     return x, y, z
pet_bar = load('data_mod.npy')
num_vowels = 10
sdwidth = 1
npts = 20
cov_mat = zeros((num_vowels,3,3))
means_mat = zeros((num_vowels,3))
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
colors = ['g','b','r','c','m','y','k','0.5']
for i in range(10):
     #change below to use different parts of the dataset
     indices = intersect1d(where( pet_bar[:,0] == 1)[0], where( pet_bar[:,2] == i+1)[0])
     # determines whether take all or >0 just takes unanimously heard correctly
     indices = intersect1d(indices, where(pet_bar[:,3] > 0.5)[0])
     pet_bar_anal = pet_bar[indices,-3:]
     cov_mat[i] = cov(pet_bar_anal, rowvar=False)
     means_mat[i] = mean(pet_bar_anal, axis=0)
     x, y, z = sphere(1, (npts,npts))
     ap = vstack((x.flatten(),y.flatten(),z.flatten()))
     d, v = eig(cov_mat[i])
     n = dot(v, (sdwidth*sqrt(d))*eye(3,3))
     out = dot(n,ap)
     bp = out + tile(means_mat[i], (npts**2,1)).T
     xp = reshape(bp[0], x.shape)
     yp = reshape(bp[1], x.shape)
     zp = reshape(bp[2], x.shape)
     ax.plot_wireframe(array(xp),array(yp),array(zp), rstride=2, cstride=2, color=colors[i%len(colors)])
ax.set_xlim3d((0,ax.get_xlim3d()[1]))
ax.set_ylim3d((0,ax.get_ylim3d()[1]))
ax.set_zlim3d((0,ax.get_zlim3d()[1]))
plt.draw()
plt.show()

3次元ガウス分布

4

1 に答える 1

1

私は基本的に、以下の完全なコードを使用して許容範囲 (ガウス球内に含まれるおおよそのポイント数) を設定し、その許容範囲での sdwidth とアルファ値を決定しました。次に、以下を使用して multivariate_normal 関数から各サンプルのアルファ値を決定できます。

temp_a = dot(dot((points-mean).T,inv(cov)),points-mean)

temp_a が決定されたアルファ値よりも小さい場合、ランダム サンプルは範囲内に収まり、受け入れられます。技術的には、ここで行われるサンプリングには受け入れ/拒否の方法を使用する必要がありますが、上記ではこれを無視しました。ほとんどの数学はここで見ることができます: http://en.wikipedia.org/wiki/Multivariate_normal_distribution

完全なコード:

gauss_toler = 0.3
value = chi2.ppf(gauss_toler,3)
lb = 1; ub = 5; runpts = 10000;
if sdwidth == None:
     sstore = -1
     sdb = 100
     alpha = 0
     i = 0
     if type(which_people) == int:
         indices = intersect1d(where( pet_bar[:,0] == 1)[0], where( pet_bar[:,2] == i+1)[0])
     else:
         indices = intersect1d(where( pet_bar[:,0] > 0)[0], where( pet_bar[:,2] == i+1)[0])
     # determines whether take all or >0 just takes unanimously heard correctly
     if unanimous_only == True:
         indices = intersect1d(indices, where(pet_bar[:,3] > 0.5)[0])
     pet_bar_anal = pet_bar[indices,-3:]
     cov_mat[i] = cov(pet_bar_anal, rowvar=False)
     means_mat[i] = mean(pet_bar_anal, axis=0)
     x, y, z = sphere(1, (npts,npts))
     ap = vstack((x.flatten(),y.flatten(),z.flatten()))
     d, v = eig(cov_mat[i])
     for sdwidth in linspace(lb,ub,runpts):
          n = dot(v, (sdwidth*sqrt(d))*eye(3,3))
          out = dot(n,ap)
          bp = out + tile(means_mat[i], (npts**2,1)).T
          xp = points_mat[i,0] = reshape(bp[0], x.shape)
          yp = points_mat[i,1] = reshape(bp[1], x.shape)
          zp = points_mat[i,2] = reshape(bp[2], x.shape)
          poi = array((points_mat[i,0,0,0], points_mat[i,1,0,0], points_mat[i,2,0,0]))
          temp_cov = cov_mat[i]
          temp_me = means_mat[i]
          a = dot(dot((poi-temp_me).T,inv(temp_cov)),poi-temp_me)
          if abs(a-value) < sdb:
             sstore = sdwidth
             alpha = a
             sdb = abs(a-value)
     sdwidth = sstore
于 2013-05-05T05:19:39.207 に答える