python - 如何将等高线图限制为二维高斯的 3sigma 变化？

Question

我对绘制二维高斯分布的等高线图有点好奇。在我的例子中，对于一组给定的 2D 点，我将它们聚集到不同的网格单元中，并计算每个单元的协方差矩阵，并绘制每个单元的高斯分布。当我绘图时，我不想要单元格的整个轮廓，但分布限制在数据点的 3 sigma 内。无论如何它可以做到吗？

我的代码如下：

import numpy as np
import matplotlib.pyplot as plt

def createCells():
    partition = 4
    coords = [np.linspace(-1.0 , 1.0, num = partition + 1) for i in range(2)]
    x, y =  np.meshgrid(*coords)
    return x, y
    
def probab(mean, covMat, lPoints):
    lPoints = lPoints[..., np.newaxis] if lPoints.ndim == 2 else lPoints  ## Create vectorized values for the x, y
    if np.linalg.det(covMat) > 0:
        factor1 = (2*np.pi)*(np.linalg.det(covMat)**(-1/2))
        factor2 = np.exp((-1/2)*np.einsum('ijk,jl,ilk->ik', lPoints - mean, np.linalg.inv(covMat), lPoints - mean))
        return factor1*factor2
    
if __name__ == '__main__':
    points = np.array([[-0.35, -0.15], [0.1, 0.1], [-0.1, 0.1], [0.05, 0.05],[0.25, 0.05], [0.1, 0.15], [0.1, 0.2], [-0.2, -0.2], [-0.25, 0.25], [0.45, 0.45], [0.75, 0.75], [0.6, 0.6], [0.55, 0.55], [0.7, 0.7], [0.68, 0.73]])
    x1, y1 = createCells()
    x = x1[0]
    y = y1[:,0]
    lP = np.array([])
    numberOftimes = 0
    for i in range(len(x) - 1):
        for j in range(len(y) - 1):
            count = 0
            meanX = 0.0
            meanY = 0.0
            localPoints = []
            covMat1 = np.array([])
            covMat2 = np.array([])
            for point in points:
                inbetween_x = x[i] <= point[0] <= x[i + 1]
                inbetween_y = y[j] <= point[1] <= y[j + 1]
                if inbetween_x and inbetween_y:
                    count += 1
                    meanX += point[0]
                    meanY += point[1]
                    localPoints.append([point[0], point[1]])
            if count >= 2:
                numberOftimes += 1
                #print(f"The local points are {localPoints}")
                localPoints = np.array(localPoints)
                meanX /= count
                meanY /= count
                meanXY = np.array([meanX, meanY])
                #print(meanXY.shape)
                #print(localPoints.shape)
                lP = localPoints - meanXY
                for k in range(count):
                    lPtranspose = (np.array([lP[k]])).T
                    lPCurrent = (np.array([lP[k]]))
                    if len(covMat1) > 0:
                        covMat1 += lPtranspose.dot(lPCurrent)
                    else:
                        covMat1 = lPtranspose*lP[k]
                        covMat1 /= count
                lPoints = localPoints[..., np.newaxis] if lP.ndim == 2 else lP  ## Create vectorized values for the x, y
                meanXY1 = localPoints.mean(0)
                meanXY2 = lPoints.mean(0)
                covMat3 = np.einsum('ijk, ikj->jk', lPoints - meanXY2, lPoints - meanXY2) / lPoints[0] - 1
                #yamlStatus = self.savingYaml(i, j, meanXY, covMat3) ## To store the cell parameters in a yaml file (for now its just out of scope for the question)
                if np.linalg.det(covMat3) > 0:   #compute the probability only if the det is not 0
                    Xx = np.linspace(x[i], x[i + 1], 1000)
                    Yy = np.linspace(y[i], y[i + 1], 1000)
                    Xx,Yy = np.meshgrid(Xx, Yy)
                    lPoints = np.vstack((Xx.flatten(), Yy.flatten())).T
                    pos = np.empty(Xx.shape + (2,))
                    pos[:, :, 0] = Xx
                    pos[:, :, 1] = Yy
                    z2 = probab(meanXY2, covMat3, lPoints)
                    summed = np.sum(z2)
                    z2 = z2.reshape(1000, 1000)
                    cs = plt.contourf(Xx, Yy, z2)#, cmap=cm.viridis)
                    plt.clabel(cs)
                localPoints = []
                #print(f"The number of times count is greater than 1 is {numberOftimes}")
        plt.plot(x1, y1, marker='.', linestyle='none', markersize=20)
        plt.plot(points[:,0 ], points[:, 1], marker='.', linestyle='none', markersize=10)
        plt.grid(linewidth=0.5)#abs(x1[0][0]-y1[0][0]))
        plt.show()

标签： pythonmatplotlibcontournormal-distribution