import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from scipy.stats import multivariate_normal # 本问题考虑的是高斯混合模型，所以导入多元高斯分布multivariate_normal
from sklearn import datasets # iris数据集加载
import math
from mpl_toolkits.mplot3d import Axes3D # 有时没这句会报错
from sklearn.svm import SVC
from sklearn.metrics import accuracy_score
from scipy.interpolate import make_interp_spline
rgb = [[255, 238, 255], # red
[228, 228, 228], # blue
[219, 238, 244]] # black
rgb = np.array(rgb ) /255.
cmap_light = ListedColormap(rgb)
iris = datasets.load_iris()
# print(iris.data) # 展示数据
def drawScatterData(iris):
x = iris.data
plt.scatter(x[0:50, 0], x[0:50, 1], color='red', marker='o', label='iris-setosa')
plt.scatter(x[50:100, 0], x[50:100, 1], color='blue', marker='*', label='iris-versicolor')
# plt.scatter(x[100:150, 0], x[100:150, 1], color='green', marker='3', label='iris-virginica')
plt.xlabel('Sepal length [cm] ')
plt.ylabel('petal length [cm] ')
plt.show()
# 统计概率
for i in range(150):
x[i, 0] = int(x[i, 0] * 10)
x[i, 1] = int(x[i, 1] * 10)
x[i, 2] = int(x[i, 2] * 10)
z = np.zeros(shape=(80, 45))
for i in range(50):
for j in range(50):
if x[i, 0] == x[j, 0] and x[i, 1] == x[j, 1]:
t1 = x[i, 0]
t1 = int(t1)
t2 = x[i, 1]
t2 = int(t2)
z[t1, t2] = z[t1, t2] + 1
for i in range(80):
for j in range(45):
z[i, j] = np.sqrt(z[i, j])
z2 = np.zeros(shape=(80, 45))
for i in range(50, 100):
for j in range(50, 100):
if x[i, 0] == x[j, 0] and x[i, 1] == x[j, 1]:
t1 = x[i, 0]
t1 = int(t1)
t2 = x[i, 1]
t2 = int(t2)
z2[t1, t2] = z2[t1, t2] + 1
for i in range(80):
for j in range(45):
z[i, j] = np.sqrt(z[i, j])
# 统计结束
x2 = np.arange(0, 8, step=0.1)
y2 = np.arange(0, 4.5, step=0.1)
xx, yy = np.meshgrid(x2, y2)
xpos, ypos = xx.ravel(), yy.ravel()
zpos = np.zeros_like(xpos)
dx = dy = 0.08
x3 = x[:, 0].ravel
y3 = x[:, 1].ravel
dz = z[:, :].ravel()
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')#三维坐标轴
ax.bar3d(xpos, ypos, zpos, dx, dy, dz, color='green', shade=True)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.show()
dz2 = z2[:, :].ravel()
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d') # 三维坐标轴
ax.bar3d(xpos, ypos, zpos, dx, dy, dz2, color='blue', shade=True)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.show()
alpha_1 = alpha_2 = 0.5
P_C1_theta = alpha_1
P_C2_theta = alpha_2
u_1 = np.zeros(shape=(2, 1))
u_1[0] = 5
u_1[1] = 3.5
u_2 = np.zeros(shape=(2, 1))
u_2[0] = 6
u_2[1] = 2.8
sigma_1 = np.zeros(shape=(2, 2))
sigma_1[0, 0] = 1
sigma_1[0, 1] = 0
sigma_1[1, 0] = 0
sigma_1[1, 1] = 1
sigma_2 = np.zeros(shape=(2, 2))
sigma_2 = sigma_1
# print(sigma_2.T) #print（转置）
# print(np.linalg.inv(sigma_2)) #print（矩阵求逆）
# print(np.linalg.norm(sigma_2)) #求模
t = iris.data
# print(np.sort(t, axis=0)) # 升序排列
# t = np.sort(t, axis=0)
x = np.zeros(shape=(2, 100))
x[0, 0:100] = t[0:100, 0].T
x[1, 0:100] = t[0:100, 1].T
# f_x_y_x_C1_theta = np.zeros(shape=(2, 150))
f_x_y_x_C1_theta = np.zeros(shape=(80, 60))
f_x_y_x_C2_theta = np.zeros(shape=(80, 60))
k_1_T = (x - u_1).T
k_1 = (x - u_1)
k_2_T = (x - u_2).T
k_2 = (x - u_2)
# 初始化均匀数据点
f_x_y_x_C1_theta_original = np.zeros(shape=(100, 100))
f_x_y_x_C2_theta_original = np.zeros(shape=(100, 100))
x_original = np.arange(0, 10, step=0.1)
y_original = np.arange(0, 10, step=0.1)
x_y_original = np.zeros(shape=(2, 100))
x_y_original[0, :] = x_original
x_y_original[1, :] = y_original
x_y_original_u_1_T = (x_y_original - u_1).T
x_y_original_u_1 = (x_y_original - u_1)
x_y_original_u_2_T = (x_y_original - u_2).T
x_y_original_u_2 = (x_y_original - u_2)
# 初始化数据点结束
for i in range(100):
# 真实数据带入计算
k_1_i_T = np.zeros(shape=(1, 2))
k_1_i_T[0, 0] = k_1_T[i, 0]
k_1_i_T[0, 1] = k_1_T[i, 1]
k_2_i_T = np.zeros(shape=(1, 2))
k_2_i_T[0, 0] = k_2_T[i, 0]
k_2_i_T[0, 1] = k_2_T[i, 1]
k_1_i = np.zeros(shape=(2, 1))
k_1_i[0, 0] = k_1[0, i]
k_1_i[1, 0] = k_1[1, i]
k_2_i = np.zeros(shape=(2, 1))
k_2_i[0, 0] = k_2[0, i]
k_2_i[1, 0] = k_2[1, i]
C_1 = (np.linalg.inv(sigma_1))
C_2 = (np.linalg.inv(sigma_2))
number_x = int(t[i, 0] * 10)
number_y = int(t[i, 1] * 10)
f_x_y_x_C1_theta[number_x, number_y] = ((np.exp(-1 / 2 * np.mat(k_1_i_T) * np.mat(C_1) * np.mat(k_1_i)))) / (
np.sqrt(2 * np.pi * 2 * np.pi) * np.linalg.norm(sigma_1))
f_x_y_x_C2_theta[number_x, number_y] = ((np.exp(-1 / 2 * np.mat(k_2_i_T) * np.mat(C_2) * np.mat(k_2_i)))) / (
np.sqrt(2 * np.pi * 2 * np.pi) * np.linalg.norm(sigma_2))
# 真实数据带入计算结束
# 初始数据带入计算
x_y_original_u_1_T_i = np.zeros(shape=(1, 2))
x_y_original_u_1_T_i[0, 0] = x_y_original_u_1_T[i, 0]
x_y_original_u_1_T_i[0, 1] = x_y_original_u_1_T[i, 1]
x_y_original_u_2_T_i = np.zeros(shape=(1, 2))
x_y_original_u_2_T_i[0, 0] = x_y_original_u_2_T[i, 0]
x_y_original_u_2_T_i[0, 1] = x_y_original_u_2_T[i, 1]
x_y_original_u_1_i = np.zeros(shape=(2, 1))
x_y_original_u_1_i[0, 0] = x_y_original_u_1[0, i]
x_y_original_u_1_i[1, 0] = x_y_original_u_1[1, i]
x_y_original_u_2_i = np.zeros(shape=(2, 1))
x_y_original_u_2_i[0, 0] = x_y_original_u_2[0, i]
x_y_original_u_2_i[1, 0] = x_y_original_u_2[1, i]
C_1 = (np.linalg.inv(sigma_1))
C_2 = (np.linalg.inv(sigma_2))
number_x = int(x_y_original[0, i] * 10)
number_y = int(x_y_original[1, i] * 10)
f_x_y_x_C1_theta_original[number_x, number_y] = ((np.exp(-1 / 2 * np.mat(x_y_original_u_1_T_i) * np.mat(C_1) * np.mat(x_y_original_u_1_i)))) / (
np.sqrt(2 * np.pi * 2 * np.pi * np.linalg.norm(sigma_1)))
f_x_y_x_C2_theta_original[number_x, number_y] = ((np.exp(-1 / 2 * np.mat(x_y_original_u_2_T_i) * np.mat(C_2) * np.mat(x_y_original_u_2_i)))) / (
np.sqrt(2 * np.pi * 2 * np.pi * np.linalg.norm(sigma_2)))
# 初始数据带入计算结束
for i in range(100):
for j in range(100):
x_y_original_u_1_T_i = np.zeros(shape=(1, 2))
x_y_original_u_1_T_i[0, 0] = x_y_original_u_1_T[i, 0]
x_y_original_u_1_T_i[0, 1] = x_y_original_u_1_T[j, 1]
x_y_original_u_2_T_i = np.zeros(shape=(1, 2))
x_y_original_u_2_T_i[0, 0] = x_y_original_u_2_T[i, 0]
x_y_original_u_2_T_i[0, 1] = x_y_original_u_2_T[j, 1]
x_y_original_u_1_i = np.zeros(shape=(2, 1))
x_y_original_u_1_i[0, 0] = x_y_original_u_1[0, i]
x_y_original_u_1_i[1, 0] = x_y_original_u_1[1, j]
x_y_original_u_2_i = np.zeros(shape=(2, 1))
x_y_original_u_2_i[0, 0] = x_y_original_u_2[0, i]
x_y_original_u_2_i[1, 0] = x_y_original_u_2[1, j]
C_1 = (np.linalg.inv(sigma_1))
C_2 = (np.linalg.inv(sigma_2))
number_x = int(x_y_original[0, i] * 10)
number_y = int(x_y_original[1, j] * 10)
f_x_y_x_C1_theta_original[number_x, number_y] = ((
np.exp(-1 / 2 * np.mat(x_y_original_u_1_T_i) * np.mat(C_1) * np.mat(x_y_original_u_1_i)))) / (
np.sqrt(2 * np.pi * 2 * np.pi * np.linalg.norm(sigma_1)))
f_x_y_x
