從高斯正態多元分佈中的某些概率繪圖numpy

Question

我有下面列出的輸出的代碼和圖片，但我想從已繪製的特定標準偏差內的這些球體中隨機抽樣。變量sdwith用於在wiremesh輸出的代碼中設置它。 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和α值。然後可以使用以下從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