• wulnch
    了解作者
  • Python
    开发工具
  • 1KB
    文件大小
  • rar
    文件格式
  • 0
    收藏次数
  • 1 积分
    下载积分
  • 0
    下载次数
  • 2021-04-14 08:45
    上传日期
a location-dispersion ellipsoid ball
test.rar
  • test.py
    2.3KB
内容介绍
# 吴凌晨 import xlrd import numpy as np import matplotlib.pyplot as plt data = xlrd.open_workbook(r'C:\\Users\\wulnc\\Desktop\\FRM\\price_data.xlsx') table = data.sheets()[3] ICBC_price = table.col_values(0) SANY_price = table.col_values(1) POLY_price = table.col_values(2) ICBC_aveprice = np.average(ICBC_price) SANY_aveprice = np.average(SANY_price) POLY_aveprice = np.average(POLY_price) col = np.array([ICBC_price, SANY_price, POLY_price]) ave = np.array([ICBC_aveprice, SANY_aveprice, POLY_aveprice]) # print(ave) cov = np.cov(col) cov_inverse = np.linalg.inv(cov) eigenvalue, eigenvector = np.linalg.eig(cov) # 特征值、特征向量 # print(cov_inverse) eigenvalue = np.sqrt(sorted(eigenvalue, reverse=True)) # eigenvalue = np.sqrt(eigenvalue) eigenvector = np.array([eigenvector[0], eigenvector[2], eigenvector[1]]) eigenvalue_matrix = np.array([[eigenvalue[0], 0, 0], [0, eigenvalue[1], 0], [0, 0, eigenvalue[2]]]) # 放入对角矩阵 # print(eigenvalue) q_lst = [] for i in col.T: a = np.dot((i-ave), cov_inverse) q_lst.append(np.dot(a, (i-ave))) q = np.average(q_lst) # print(q) # 椭球绘图 def plot(q,color): ax = plt.subplot(111, projection='3d') theta = np.linspace(0, 2 * np.pi, 100) zl = (np.linspace(0, np.pi, 15)) sin = np.sin(zl) zcos = np.cos(zl) for i in range(15): z = [(zcos[i] * q) for k in range(100)] x = sin[i] * np.sin(theta) * q y = sin[i] * np.cos(theta) * q y_coordinate = np.array([x, y, z]) # 单位球的坐标 # 坐标变换 transform = np.dot(eigenvector, eigenvalue_matrix) x_coordinate = np.dot(transform, y_coordinate) x_1 = x_coordinate[0]+ave[0] x_2 = x_coordinate[1]+ave[1] x_3 = x_coordinate[2]+ave[2] if i == 0: ax.plot(x_1, x_2, x_3, label='parametric curve') else: ax.plot(x_1, x_2, x_3,c=color) x, y, z = col[0], col[1], col[2] ax = plt.subplot(111, projection='3d') ax.scatter(x, y, z, s=5, c='g') ax.set_xlim3d((-0.08, 0.08)) ax.set_ylim3d((-0.08, 0.08)) ax.set_zlim3d((-0.08, 0.08)) ax.set_zlabel('ICBC') # 坐标轴 ax.set_ylabel('SANY') ax.set_xlabel('POLY') # plot(q,'b') plot(q*0.95,'r') plot(q*1.05,'c') plot(q*1.1,'y') plot(q*0.9,'k') plt.show()
评论
    相关推荐