如何计算椭圆高斯分布的角度 [英] How to calculate the angle of ellipse Gaussian distribution

问题描述

我使用以下Python代码来计算矩量法的高斯式分布基础的中心和大小.但是，我无法编写代码来计算高斯角.

I make following Python Code to calculate center and size of Gaussian-like distribution basis of moment method. But, I can't make the code to calculate the angle of gaussian.

请看图片.

第一张图片是原始数据.

First Picture is original data.

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

Second picture is reconstruct data from the result of moment method.

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

But, second picture is insufficient reconstruction. Because, original data is inclined distribution. I have to, I think, calculate the angle of axis for Gaussian-like distribution.

假设原始分布足够像高斯分布.

To assume that the original distribution is sufficiently Gaussian-like distribution.

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

推荐答案

正如Paul Panzer所说，该方法的缺点是您寻找权重"和角度"而不是协方差矩阵.协方差矩阵非常适合您的方法:只需再计算一个矩，混合xy即可.

As Paul Panzer said, the flaw of your approach is that you look for "weight" and "angle" instead of covariance matrix. The covariance matrix fits perfectly in your approach: just compute one more moment, mixed 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，现在可以通过以下方式重建它:

for reconstruction purpose. Still using your make_gauss function to create the original PDF, this is how it now gets reconstructed:

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

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

That's it; reconstruction works fine now.

颜色栏上的单位不同，因为您的 make_gauss 公式无法规范PDF.

Units on the color bar are not the same, because your make_gauss formula does not normalize the PDF.

这篇关于如何计算椭圆高斯分布的角度的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！