Блог

Мне было немного любопытно построить контурные графики для 2-мерного распределения Гаусса. В моем случае для заданного набора 2D-точек я группирую их в разные ячейки сетки и вычисляю матрицу ковариации для каждой ячейки и строю распределение Гаусса для каждой ячейки. Когда я строю график, мне нужен не весь контур ячейки, а распределение, ограниченное в пределах 3 сигм точек данных. Есть ли вообще возможность это сделать ?

Мой код выглядит следующим образом:

 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()