python - 从 numpy 中高斯正态多元分布的某些概率中得出

Question

我有下面列出的输出的代码和图片，但我想在已绘制的特定标准差范围内从这些球体中抽取随机样本。变量 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()

Answer 1

我基本上使用下面的完整代码来设置容差(高斯球内包含的大约点数)以确定该容差下的 sdwidth 和 alpha 值。然后可以使用以下方法从 multivariate_normal 函数确定每个样本的 alpha 值:

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

如果 temp_a 小于确定的 alpha 值，则随机样本落入球体内并被接受。值得注意的是，从技术上讲，接受/拒绝方法应该用于此处完成的采样，但我在上面忽略了这一点。大多数数学可以在这里查看: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