Как рассчитать угол эллипса гауссовского распределения

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

Как сказал Пол Панцер, недостатком вашего подхода является то, что вы ищете "вес" и "угол" вместо ковариационной матрицы. Ковариационная матрица идеально подходит для вашего подхода: просто вычислите еще один момент, смешанный 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.