使用matplotlib绘制比例缩放和旋转的双变量分布 [英] Plot scaled and rotated bivariate distribution using matplotlib

查看:142
本文介绍了使用matplotlib绘制比例缩放和旋转的双变量分布的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

我正在尝试绘制绘图一个双变量 高斯使用 matplotlib 发行版。我想使用两个散点点(组A)(组B)的 xy 坐标来执行此操作。

I am trying to plot a bivariate gaussian distribution using matplotlib. I want to do this using the xy coordinates of two scatter points (Group A), (Group B).

我想通过调整 COV 分布 > 矩阵来说明每个组速度以及它们与另一个 xy 坐标用作参考点。

I want to adjust the distribution by adjusting the COV matrix to account for each Groups velocity and their distance to an additional xy coordinate used as a reference point.

我已计算出每个组 xy 坐标与参考点的距离。距离以半径表示,标记为 [GrA_Rad] [GrB_Rad]

I've calculated the distance of each groups xy coordinate to that of the reference point. The distance is expressed as a radius, labelled [GrA_Rad],[GrB_Rad].

所以它们离参考点越远,半径越大。我还计算了标记为 [GrA_Vel] [GrB_Vel] 速度 $ c>。每个组的方向表示为方向。这被标记为 [GrA_Rotation] [GrB_Rotation]

So the further they are away from the reference point the greater the radius. I've also calculated velocity labelled [GrA_Vel],[GrB_Vel]. The direction of each group is expressed as the orientation. This is labelled [GrA_Rotation],[GrB_Rotation]

有关我希望如何针对速度和距离调整分布的问题(半径)

Question on how I want the distribution to be adjusted for velocity and distance (radius):

我希望使用 SVD 。具体来说,如果每个散点的旋转角度为 ,这提供了方向速度可用于描述缩放比例 矩阵 [GrA_Scaling] [GrB_Scaling] 。因此,此缩放比例 矩阵可用于扩展 radius 沿 x方向并沿 y方向 radius $ c>。这表示 COV 矩阵

I'm hoping to use SVD. Specifically, if I have the rotation angle of each scatter, this provides the direction. The velocity can be used to describe a scaling matrix [GrA_Scaling],[GrB_Scaling]. So this scaling matrix can be used to expand the radius in the x-direction and contract the radius in the y-direction. This expresses the COV matrix.

最后,找到分布 平均值的值是将组位置 (x,y)转换为速度

Finally, the distribution mean value is found by translating the groups location (x,y) by half the velocity.

简单地输入半径应用于每个组的散点点。 COV矩阵通过半径速度进行调整。因此,使用缩放比例 矩阵来扩展 radius x方向并沿 y方向收缩。 方向是从旋转 角度测量的。然后通过翻译组位置(x,y)分布 平均值的值是速度 的一半。

Put simply: the radius is applied to each group's scatter point. The COV matrix is adjusted by the radius and velocity. So using the scaling matrix to expand the radius in x-direction and contract in y-direction. The direction is measured from the rotation angle. Then determine the distribution mean value by translating the groups location (x,y) by half the velocity.

下面是 df 这些变量

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.animation as animation

d = ({
    'Time' : [1,2,3,4,5,6,7,8],       
    'GrA_X' : [10,12,17,16,16,14,12,8],                 
    'GrA_Y' : [10,12,13,7,6,7,8,8], 
    'GrB_X' : [5,8,13,16,19,15,13,5],                 
    'GrB_Y' : [6,15,12,7,8,9,10,8],   
    'Reference_X' : [6,8,14,18,13,11,16,15],                 
    'Reference_Y' : [10,12,8,12,15,12,10,8],                  
    'GrA_Rad' : [8.3,8.25,8.2,8,8.15,8.15,8.2,8.3],  
    'GrB_Rad' : [8.3,8.25,8.3,8.4,8.6,8.4,8.3,8.65],               
    'GrA_Vel' : [0,2.8,5.1,6.1,1.0,2.2,2.2,4.0],
    'GrB_Vel' : [0,9.5,5.8,5.8,3.16,4.12,2.2,8.2],               
    'GrA_Scaling' : [0,0.22,0.39,0.47,0.07,0.17,0.17,0.31],
    'GrB_Scaling' : [0,0.53,0.2,0.2,0.06,0.1,0.03,0.4],                   
    'GrA_Rotation' : [0,45,23.2,-26.56,-33.69,-36.86,-45,-135], 
    'GrB_Rotation' : [0,71.6,36.87,5.2,8.13,16.70,26.57,90],                       
     })

df = pd.DataFrame(data = d)

我已经制作了动画 每个 xy坐标

I've made an animated plot of each xy coordinate.

GrA_X = [10,12,17,16,16,14,12,8]
GrA_Y = [10,12,13,7,6,7,8,8]

GrB_X = [5,8,13,16,19,15,13,5]                 
GrB_Y = [6,15,12,10,8,9,10,8]

Item_X = [6,8,14,18,13,11,16,15]  
Item_Y = [10,12,8,12,15,12,10,8]

scatter_GrA = ax.scatter(GrA_X, GrA_Y) 
scatter_GrB = ax.scatter(GrB_X, GrB_Y) 
scatter_Item = ax.scatter(Item_X, Item_Y) 

def animate(i) :
    scatter_GrA.set_offsets([[GrA_X[0+i], GrA_Y[0+i]]])
    scatter_GrB.set_offsets([[GrB_X[0+i], GrB_Y[0+i]]])
    scatter_Item.set_offsets([[Item_X[0+i], Item_Y[0+i]]])    

ani = animation.FuncAnimation(fig, animate, np.arange(0,9),
                              interval = 1000, blit = False)


推荐答案

更新



问题已更新,并且变得更加清晰。我已更新代码以使其匹配。这是最新的输出:

Update

The question has been updated, and has gotten somewhat clearer. I've updated my code to match. Here's the latest output:

以下是用于生成上述情节的代码:

Here's the code that was used to produce the above plot:

dfake = ({    
    'GrA_X' : [15,15],                 
    'GrA_Y' : [15,15], 
    'Reference_X' : [15,3],                 
    'Reference_Y' : [15,15],                  
    'GrA_Rad' : [15,25],                 
    'GrA_Vel' : [0,10],
    'GrA_Scaling' : [0,0.5],
    'GrA_Rotation' : [0,45]                     
})

dffake = pd.DataFrame(dfake)
fig,axs = plt.subplots(1, 2, figsize=(16,8))
fig.subplots_adjust(0,0,1,1)
plotone(dffake, 'A', 0, xlim=(0,30), ylim=(0,30), fig=fig, ax=axs[0])
plotone(dffake, 'A', 1, xlim=(0,30), ylim=(0,30), fig=fig, ax=axs[1])
plt.show()

和我使用的 plotone 函数的完整实现在下面的代码块中。如果您只是想了解用于生成和转换2D高斯PDF的数学方法,请查看 mvpdf 函数(以及 rot getcov 函数所依赖):

and the complete implementation of the plotone function that I used is in the code block below. If you just want to know about the math used to generate and transform the 2D gaussian PDF, check out the mvpdf function (and the rot and getcov functions it depends on):

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy.stats as sts

def rot(theta):
    theta = np.deg2rad(theta)
    return np.array([
        [np.cos(theta), -np.sin(theta)],
        [np.sin(theta), np.cos(theta)]
    ])

def getcov(radius=1, scale=1, theta=0):
    cov = np.array([
        [radius*(scale + 1), 0],
        [0, radius/(scale + 1)]
    ])

    r = rot(theta)
    return r @ cov @ r.T

def mvpdf(x, y, xlim, ylim, radius=1, velocity=0, scale=0, theta=0):
    """Creates a grid of data that represents the PDF of a multivariate gaussian.

    x, y: The center of the returned PDF
    (xy)lim: The extent of the returned PDF
    radius: The PDF will be dilated by this factor
    scale: The PDF be stretched by a factor of (scale + 1) in the x direction, and squashed by a factor of 1/(scale + 1) in the y direction
    theta: The PDF will be rotated by this many degrees

    returns: X, Y, PDF. X and Y hold the coordinates of the PDF.
    """
    # create the coordinate grids
    X,Y = np.meshgrid(np.linspace(*xlim), np.linspace(*ylim))

    # stack them into the format expected by the multivariate pdf
    XY = np.stack([X, Y], 2)

    # displace xy by half the velocity
    x,y = rot(theta) @ (velocity/2, 0) + (x, y)

    # get the covariance matrix with the appropriate transforms
    cov = getcov(radius=radius, scale=scale, theta=theta)

    # generate the data grid that represents the PDF
    PDF = sts.multivariate_normal([x, y], cov).pdf(XY)

    return X, Y, PDF

def plotmv(x, y, xlim=None, ylim=None, radius=1, velocity=0, scale=0, theta=0, xref=None, yref=None, fig=None, ax=None):
    """Plot an xy point with an appropriately tranformed 2D gaussian around it.
    Also plots other related data like the reference point.
    """
    if xlim is None: xlim = (x - 5, x + 5)
    if ylim is None: ylim = (y - 5, y + 5)

    if fig is None:
        fig = plt.figure(figsize=(8,8))
        ax = fig.gca()
    elif ax is None:
        ax = fig.gca()

    # plot the xy point
    ax.plot(x, y, '.', c='C0', ms=20)

    if not (xref is None or yref is None):
        # plot the reference point, if supplied
        ax.plot(xref, yref, '.', c='w', ms=12)

    # plot the arrow leading from the xy point
    if velocity > 0:
        ax.arrow(x, y, *rot(theta) @ (velocity, 0), 
                 width=.4, length_includes_head=True, ec='C0', fc='C0')

    # fetch the PDF of the 2D gaussian
    X, Y, PDF = mvpdf(x, y, xlim=xlim, ylim=ylim, radius=radius, velocity=velocity, scale=scale, theta=theta)

    # normalize PDF by shifting and scaling, so that the smallest value is 0 and the largest is 1
    normPDF = PDF - PDF.min()
    normPDF = normPDF/normPDF.max()

    # plot and label the contour lines of the 2D gaussian
    cs = ax.contour(X, Y, normPDF, levels=6, colors='w', alpha=.5)
    ax.clabel(cs, fmt='%.3f', fontsize=12)

    # plot the filled contours of the 2D gaussian. Set levels high for smooth contours
    cfs = ax.contourf(X, Y, normPDF, levels=50, cmap='viridis', vmin=-.9, vmax=1)

    # create the colorbar and ensure that it goes from 0 -> 1
    cbar = fig.colorbar(cfs, ax=ax)
    cbar.set_ticks([0, .2, .4, .6, .8, 1])

    # add some labels
    ax.grid()
    ax.set_xlabel('X distance (M)')
    ax.set_ylabel('Y distance (M)')

    # ensure that x vs y scaling doesn't disrupt the transforms applied to the 2D gaussian
    ax.set_aspect('equal', 'box')

    return fig, ax

def fetchone(df, l, i, **kwargs):
    """Fetch all the needed data for one xy point
    """
    keytups = (
        ('x', 'Gr%s_X'%l),
        ('y', 'Gr%s_Y'%l),
        ('radius', 'Gr%s_Rad'%l),
        ('velocity', 'Gr%s_Vel'%l),
        ('scale', 'Gr%s_Scaling'%l),
        ('theta', 'Gr%s_Rotation'%l),
        ('xref', 'Reference_X'),
        ('yref', 'Reference_Y')
    )

    ret = {k:df.loc[i, l] for k,l in keytups}
    # add in any overrides
    ret.update(kwargs)

    return ret

def plotone(df, l, i, xlim=None, ylim=None, fig=None, ax=None, **kwargs):
    """Plot exactly one point from the dataset
    """
    # look up all the data to plot one datapoint
    xydata = fetchone(df, l, i, **kwargs)

    # do the plot
    return plotmv(xlim=xlim, ylim=ylim, fig=fig, ax=ax, **xydata)



旧答案-2



我已经调整了答案以匹配OP发布的示例:

Old answer -2

I've adjusted my answer to match the example the OP posted:

< a href = https://i.stack.imgur.com/t39Ud.png rel = nofollow noreferrer>

以下是产生以上图像的代码:

Here's the code that produced the above image:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy.stats as sts

def rot(theta):
    theta = np.deg2rad(theta)
    return np.array([
        [np.cos(theta), -np.sin(theta)],
        [np.sin(theta), np.cos(theta)]
    ])

def getcov(radius=1, scale=1, theta=0):
    cov = np.array([
        [radius*(scale + 1), 0],
        [0, radius/(scale + 1)]
    ])

    r = rot(theta)
    return r @ cov @ r.T

def datalimits(*data, pad=.15):
    dmin,dmax = min(d.min() for d in data), max(d.max() for d in data)
    spad = pad*(dmax - dmin)
    return dmin - spad, dmax + spad

d = ({
    'Time' : [1,2,3,4,5,6,7,8],       
    'GrA_X' : [10,12,17,16,16,14,12,8],                 
    'GrA_Y' : [10,12,13,7,6,7,8,8], 
    'GrB_X' : [5,8,13,16,19,15,13,5],                 
    'GrB_Y' : [6,15,12,7,8,9,10,8],   
    'Reference_X' : [6,8,14,18,13,11,16,15],                 
    'Reference_Y' : [10,12,8,12,15,12,10,8],                  
    'GrA_Rad' : [8.3,8.25,8.2,8,8.15,8.15,8.2,8.3],  
    'GrB_Rad' : [8.3,8.25,8.3,8.4,8.6,8.4,8.3,8.65],               
    'GrA_Vel' : [0,2.8,5.1,6.1,1.0,2.2,2.2,4.0],
    'GrB_Vel' : [0,9.5,5.8,5.8,3.16,4.12,2.2,8.2],               
    'GrA_Scaling' : [0,0.22,0.39,0.47,0.07,0.17,0.17,0.31],
    'GrB_Scaling' : [0,0.53,0.2,0.2,0.06,0.1,0.03,0.4],                   
    'GrA_Rotation' : [0,45,23.2,-26.56,-33.69,-36.86,-45,-135], 
    'GrB_Rotation' : [0,71.6,36.87,5.2,8.13,16.70,26.57,90],                       
     })

df = pd.DataFrame(data=d)

limitpad = .5
clevels = 5
cflevels = 50

xmin,xmax = datalimits(df['GrA_X'], df['GrB_X'], pad=limitpad)
ymin,ymax = datalimits(df['GrA_Y'], df['GrB_Y'], pad=limitpad)

X,Y = np.meshgrid(np.linspace(xmin, xmax), np.linspace(ymin, ymax))

fig = plt.figure(figsize=(10,6))
ax = plt.gca()

Zs = []
for l,color in zip('AB', ('red', 'yellow')):
    # plot all of the points from a single group
    ax.plot(df['Gr%s_X'%l], df['Gr%s_Y'%l], '.', c=color, ms=15, label=l)

    Zrows = []
    for _,row in df.iterrows():
        x,y = row['Gr%s_X'%l], row['Gr%s_Y'%l]

        cov = getcov(radius=row['Gr%s_Rad'%l], scale=row['Gr%s_Scaling'%l], theta=row['Gr%s_Rotation'%l])
        mnorm = sts.multivariate_normal([x, y], cov)
        Z = mnorm.pdf(np.stack([X, Y], 2))
        Zrows.append(Z)

    Zs.append(np.sum(Zrows, axis=0))

# plot the reference points

# create Z from the difference of the sums of the 2D Gaussians from group A and group B
Z = Zs[0] - Zs[1]

# normalize Z by shifting and scaling, so that the smallest value is 0 and the largest is 1
normZ = Z - Z.min()
normZ = normZ/normZ.max()

# plot and label the contour lines
cs = ax.contour(X, Y, normZ, levels=clevels, colors='w', alpha=.5)
ax.clabel(cs, fmt='%2.1f', colors='w')#, fontsize=14)

# plot the filled contours. Set levels high for smooth contours
cfs = ax.contourf(X, Y, normZ, levels=cflevels, cmap='viridis', vmin=0, vmax=1)
# create the colorbar and ensure that it goes from 0 -> 1
cbar = fig.colorbar(cfs, ax=ax)
cbar.set_ticks([0, .2, .4, .6, .8, 1])


ax.set_aspect('equal', 'box')



旧答案-1



要确切地说明您要追求的目标有点困难。可以通过其协方差矩阵缩放和旋转多元高斯分布。以下是根据您的数据执行此操作的示例:

Old answer -1

It's a little hard to tell exactly what you're after. It is possible to scale and rotate a multivariate gaussian distribution via its covariance matrix. Here's an example of how to do so based on your data:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy.stats as sts

def rot(theta):
    theta = np.deg2rad(theta)
    return np.array([
        [np.cos(theta), -np.sin(theta)],
        [np.sin(theta), np.cos(theta)]
    ])

def getcov(scale, theta):
    cov = np.array([
        [1*(scale + 1), 0],
        [0, 1/(scale + 1)]
    ])

    r = rot(theta)
    return r @ cov @ r.T

d = ({
    'Time' : [1,2,3,4,5,6,7,8],       
    'GrA_X' : [10,12,17,16,16,14,12,8],                 
    'GrA_Y' : [10,12,13,7,6,7,8,8], 
    'GrB_X' : [5,8,13,16,19,15,13,5],                 
    'GrB_Y' : [6,15,12,7,8,9,10,8],   
    'Reference_X' : [6,8,14,18,13,11,16,15],                 
    'Reference_Y' : [10,12,8,12,15,12,10,8],                  
    'GrA_Rad' : [8.3,8.25,8.2,8,8.15,8.15,8.2,8.3],  
    'GrB_Rad' : [8.3,8.25,8.3,8.4,8.6,8.4,8.3,8.65],               
    'GrA_Vel' : [0,2.8,5.1,6.1,1.0,2.2,2.2,4.0],
    'GrB_Vel' : [0,9.5,5.8,5.8,3.16,4.12,2.2,8.2],               
    'GrA_Scaling' : [0,0.22,0.39,0.47,0.07,0.17,0.17,0.31],
    'GrB_Scaling' : [0,0.53,0.2,0.2,0.06,0.1,0.03,0.4],                   
    'GrA_Rotation' : [0,45,23.2,-26.56,-33.69,-36.86,-45,-135], 
    'GrB_Rotation' : [0,71.6,36.87,5.2,8.13,16.70,26.57,90],                       
     })

df = pd.DataFrame(data=d)
xmin,xmax = min(df['GrA_X'].min(), df['GrB_X'].min()), max(df['GrA_X'].max(), df['GrB_X'].max())
ymin,ymax = min(df['GrA_Y'].min(), df['GrB_Y'].min()), max(df['GrA_Y'].max(), df['GrB_Y'].max())

X,Y = np.meshgrid(
    np.linspace(xmin - (xmax - xmin)*.1, xmax + (xmax - xmin)*.1),
    np.linspace(ymin - (ymax - ymin)*.1, ymax + (ymax - ymin)*.1)
)

fig,axs = plt.subplots(df.shape[0], sharex=True, figsize=(4, 4*df.shape[0]))
fig.subplots_adjust(0,0,1,1,0,-.82)

for (_,row),ax in zip(df.iterrows(), axs):
    for c in 'AB':
        x,y = row['Gr%s_X'%c], row['Gr%s_Y'%c]

        cov = getcov(scale=row['Gr%s_Scaling'%c], theta=row['Gr%s_Rotation'%c])
        mnorm = sts.multivariate_normal([x, y], cov)
        Z = mnorm.pdf(np.stack([X, Y], 2))

        ax.contour(X, Y, Z)

        ax.plot(row['Gr%s_X'%c], row['Gr%s_Y'%c], 'x')
        ax.set_aspect('equal', 'box')

输出如下:

这篇关于使用matplotlib绘制比例缩放和旋转的双变量分布的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!

查看全文
登录 关闭
扫码关注1秒登录
发送“验证码”获取 | 15天全站免登陆