python绘制图像中心坐标二维分布曲线

在这里插入图片描述

数据和代码如下所示：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import xlrd
from scipy.stats import multivariate_normal
from mpl_toolkits.mplot3d import Axes3D
np.set_printoptions(suppress=True)

# 根据均值、标准差,求指定范围的正态分布概率值
def normfun(x, mu, sigma):
  pdf = np.exp(-((x - mu)**2)/(2*sigma**2)) / (sigma * np.sqrt(2*np.pi))
  return pdf

# 读取excel表格数据
data = xlrd.open_workbook("laserspot.xls")  # 打开excel文件
table1 = data.sheet_by_name('Sheet2')  # 通过excel里面的表名获取工作表
row1 = table1.row_values(0)  # 根据索引读取一行的数据，参数：行索引，开始列索引，结束列索引（不包含）
for i in range(len(row1)):  # 获取表头
    if row1[i] == 'U':
        oneindex1 = i
    elif row1[i] == 'V':
        oneindex2 = i

# 获取U、V三列的数据
u = table1.col_values(oneindex1, 1)
v = table1.col_values(oneindex2, 1)
u = np.array(u)
v = np.array(v)

mu1 = u.mean()
std1 = u.std()
mu2 = v.mean()
std2 = v.std()

z=np.stack((u,v))
rou = np.cov(z)[0,1]
# print(rou)

# 输出U的数据
x1 = np.arange(min(u), max(u), 0.001)
# 设定 y 轴，载入刚才的正态分布函数
y1 = normfun(x1, u.mean(), u.std())

# 输出V的数据
x2 = np.arange(min(v), max(v), 0.001)
# 设定 y 轴，载入刚才的正态分布函数
y2 = normfun(x2, v.mean(), v.std())


X, Y = np.meshgrid(x1, x2)
# 二维坐标数据
d = np.dstack([X,Y])
print(d)
# 计算二维联合高斯概率
Z = 1 / ( (2*np.pi) * std1 *std2 * np.sqrt(1-rou**2)) * np.exp((-( ((X - mu1)**2 / std1**2) - (2*rou*(X-mu1)*(Y-mu2)/std1*std2) + ((Y - mu2)**2 / std2**2) ))/(2*(1- rou**2)))

'''绘制直方分布图'''
fig1 = plt.figure(figsize=(11, 10))
ax = fig1.add_subplot(2, 1, 1)
ax.plot(x1, y1)
ax.hist(u, bins=25, rwidth=0.8, density=True) # bins个柱状图,宽度是rwidth(0~1),=1没有缝隙
ax.set_title("U")
ax.set_xlabel('Coordinate')
ax.set_ylabel('probability')

ax2 = fig1.add_subplot(2, 1, 2)
ax2.plot(x2, y2)
ax2.hist(v, bins=25, rwidth=0.8, density=True) # bins个柱状图,宽度是rwidth(0~1),=1没有缝隙
ax2.set_title('V')
ax2.set_xlabel('Coordinate')
ax2.set_ylabel('probability')
ax2.ticklabel_format(useOffset=False, style='plain')

'''绘制二元高斯概率分布图'''
fig2 = plt.figure(figsize=(8,6))
ax = Axes3D(fig2)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='rainbow', alpha=0.8)
# ax.plot_trisurf(X, Y, Z, rstride=1, cstride=1, cmap='seismic', alpha=0.8)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.ticklabel_format(useOffset=False, style='plain')

plt.show()