python - 如何计算椭圆高斯分布的角度

Question

我编写了以下Python代码来计算矩量方法的类高斯分布基的中心和大小。但是，我无法编写计算高斯角度的代码。

请看图片。

第一张图是原始数据。

第二张图是根据矩量法的结果重建数据。

但是，第二张图片重建不充分。因为，原始数据是倾斜分布的。我认为，我必须计算类高斯分布的轴角。

假设原始分布足够类似高斯分布。

import numpy as np
import matplotlib.pyplot as plt
import json, glob
import sys, time, os
from mpl_toolkits.axes_grid1 import make_axes_locatable
from linecache import getline, clearcache
from scipy.integrate import simps
from scipy.constants import *

def integrate_simps (mesh, func):
    nx, ny = func.shape
    px, py = mesh[0][int(nx/2), :], mesh[1][:, int(ny/2)]
    val = simps( simps(func, px), py )
    return val

def normalize_integrate (mesh, func):
    return func / integrate_simps (mesh, func)

def moment (mesh, func, index):
    ix, iy = index[0], index[1]
    g_func = normalize_integrate (mesh, func)
    fxy = g_func * mesh[0]**ix * mesh[1]**iy
    val = integrate_simps (mesh, fxy)
    return val

def moment_seq (mesh, func, num):
    seq = np.empty ([num, num])
    for ix in range (num):
        for iy in range (num):
            seq[ix, iy] = moment (mesh, func, [ix, iy])
    return seq

def get_centroid (mesh, func):
    dx = moment (mesh, func, (1, 0))
    dy = moment (mesh, func, (0, 1))
    return dx, dy

def get_weight (mesh, func, dxy):
    g_mesh = [mesh[0]-dxy[0], mesh[1]-dxy[1]]
    lx = moment (g_mesh, func, (2, 0))
    ly = moment (g_mesh, func, (0, 2))
    return np.sqrt(lx), np.sqrt(ly)

def plot_contour_sub (mesh, func, loc=[0, 0], title="name", pngfile="./name"):
    sx, sy = loc
    nx, ny = func.shape
    xs, ys = mesh[0][0, 0], mesh[1][0, 0]
    dx, dy = mesh[0][0, 1] - mesh[0][0, 0], mesh[1][1, 0] - mesh[1][0, 0]
    mx, my = int ( (sy-ys)/dy ), int ( (sx-xs)/dx )
    fig, ax = plt.subplots()
    divider = make_axes_locatable(ax)
    ax.set_aspect('equal')
    ax_x = divider.append_axes("bottom", 1.0, pad=0.5, sharex=ax)
    ax_x.plot (mesh[0][mx, :], func[mx, :])
    ax_x.set_title ("y = {:.2f}".format(sy))
    ax_y = divider.append_axes("right" , 1.0, pad=0.5, sharey=ax)
    ax_y.plot (func[:, my], mesh[1][:, my])
    ax_y.set_title ("x = {:.2f}".format(sx))
    im = ax.contourf (*mesh, func, cmap="jet")
    ax.set_title (title)
    plt.colorbar (im, ax=ax, shrink=0.9)
    plt.savefig(pngfile + ".png")

def make_gauss (mesh, sxy, rxy, rot):
    x, y = mesh[0] - sxy[0], mesh[1] - sxy[1]
    px = x * np.cos(rot) - y * np.sin(rot)
    py = y * np.cos(rot) + x * np.sin(rot)
    fx = np.exp (-0.5 * (px/rxy[0])**2)
    fy = np.exp (-0.5 * (py/rxy[1])**2)
    return fx * fy

if __name__ == "__main__":
    argvs = sys.argv  
    argc = len(argvs)
    print (argvs)

    nx, ny = 500, 500
    lx, ly = 200, 150
    rx, ry = 40, 25
    sx, sy = 50, 10
    rot    = 30

    px = np.linspace (-1, 1, nx) * lx
    py = np.linspace (-1, 1, ny) * ly
    mesh = np.meshgrid (px, py)
    fxy0 = make_gauss (mesh, [sx, sy], [rx, ry], np.deg2rad(rot)) * 10
    s0xy = get_centroid (mesh, fxy0)
    w0xy = get_weight (mesh, fxy0, s0xy)

    fxy1 = make_gauss (mesh, s0xy, w0xy, np.deg2rad(0))
    s1xy = get_centroid (mesh, fxy1)
    w1xy = get_weight (mesh, fxy1, s1xy)

    print ([sx, sy], s0xy, s1xy)
    print ([rx, ry], w0xy, w1xy)

    plot_contour_sub (mesh, fxy0, loc=s0xy, title="Original", pngfile="./fxy0")
    plot_contour_sub (mesh, fxy1, loc=s1xy, title="Reconst" , pngfile="./fxy1")

Answer 1

正如 Paul Panzer 所说，你的方法的缺陷是你寻找“权重”和“角度”而不是协方差矩阵。协方差矩阵完全适合您的方法:只需再计算一个时刻，混合 xy。

函数 get_weight 应替换为

def get_covariance (mesh, func, dxy):
    g_mesh = [mesh[0]-dxy[0], mesh[1]-dxy[1]]
    Mxx = moment (g_mesh, func, (2, 0))
    Myy = moment (g_mesh, func, (0, 2))
    Mxy = moment (g_mesh, func, (1, 1))
    return np.array([[Mxx, Mxy], [Mxy, Myy]])

再添加一个导入，

from scipy.stats import multivariate_normal

用于重建目的。仍然使用 make_gauss 函数创建原始 PDF，这就是它现在重建的方式:

s0xy = get_centroid (mesh, fxy0)
w0xy = get_covariance (mesh, fxy0, s0xy)
fxy1 = multivariate_normal.pdf(np.stack(mesh, -1), mean=s0xy, cov=w0xy)

就是这样；现在重建工作正常。

颜色条上的单位不相同，因为您的 make_gauss 公式没有对 PDF 进行标准化。