楕円ガウス分布の角度を計算する方法

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

ポール装甲はあなたのアプローチの欠陥があなたの代わりに共分散行列の「量」と「アングル」を探していることである、言ったように。共分散行列は、あなたのアプローチに完全に合っています。

機能get_weightは

from scipy.stats import multivariate_normal 

復興目的のために、1つの以上のインポートを追加

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

に置き換える必要があります。それでも、元のPDFを作成するためにあなたのmake_gauss機能を使って、これは今で再構成します方法です。それだけです

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を正規化していないため、カラーバー上の

ユニットは、同じではありません。