『科學計算』可視化二元正態分佈&3D科學可視化實戰

二元正態分佈可視化本體 

因爲近來一直再看kaggle的入門書（sklearn入門手冊的感受233），感受對機器學習的理解加深了很多（實際上就只是調包能力增強了），聯想到假期在python科學計算上也算是進行了一些嘗試學習，以爲仍是須要學習一下機器學習原理的，因此從新啃起了吳恩達的cs229，上次（5月份的時候？）就是在多元高斯分佈這裏吃的癟，看不下去了，此次覺定穩紮穩打，不求速度多實踐實踐，儘可能理解數學原理，因此再次看到這部分時決定把這個分佈復現出來，吳恩達大佬用的matlab，我用的python，畫的還不錯，代碼以下，

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
from matplotlib import cm
import matplotlib as mpl

num = 200
l = np.linspace(-5,5,num)
X, Y =np.meshgrid(l, l)

u = np.array([0, 0])
o = np.array([[1, 0.5],
              [0.5, 1]])

pos = np.concatenate((np.expand_dims(X,axis=2),np.expand_dims(Y,axis=2)),axis=2)

a = (pos-u).dot(np.linalg.inv(o))
b = np.expand_dims(pos-u,axis=3)
Z = np.zeros((num,num), dtype=np.float32)
for i in range(num):
    Z[i] = [np.dot(a[i,j],b[i,j]) for j in range(num)]
Z = np.exp(Z*(-0.5))/(2*np.pi*np.linalg.det(o))
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
ax.plot_surface(X, Y, Z, rstride=5, cstride=5, alpha=0.3, cmap=cm.coolwarm)

cset = ax.contour(X,Y,Z,10,zdir='z',offset=0,cmap=cm.coolwarm)
cset = ax.contour(X, Y, Z, zdir='x', offset=-5,cmap=mpl.cm.winter)
cset = ax.contour(X, Y, Z, zdir='y', offset= 5,cmap= mpl.cm.winter)
'''
mpl.cm.rainbow
mpl.cm.winter
mpl.cm.bwr  # 藍，白，紅
cm.coolwarm
'''

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()　

實際操做中，能夠看到我在Z生成部分使用了雙層循環，我本意是使用numpy廣播機制優化掉循環，實際操做不太順利，（20，20，2）去叉乘（20，20，2，1），結果shape不是我指望的（20，20，1），而是（20，20，20，20，1），也就是說在高維叉乘時其實廣播機制不太好用，畢竟實際上兩個不一樣維度矩陣是能夠直接叉乘的（雖然對維度有要求），這一點值得注意（高維矩陣叉乘不要依賴numpy的廣播機制）。

 

參數：

u = np.array([0, 0]) 
o = np.array([[1, 0.5], 
             [0.5, 1]])

　

參數：

u = np.array([1, 1])
o = np.array([[1, 0],
              [0, 1]])

 參數：

u = np.array([1, 1])
o = 3*np.array([[1, 0],
              [0, 1]])

 

 

等高線圖添加

咱們單獨繪製一下等高線圖，

# 前面添加圖的位置修改以下，
# ax = fig.add_subplot(211,projection='3d')

ax2 = fig.add_subplot(212)
cs = ax2.contour(X,Y,Z)
ax2.clabel(cs, inline=1, fontsize=20)

高斯判別分析模型示意圖可視化

如今咱們在上面代碼的基礎上可視化吳恩達老大的下一節的圖，高斯判別分析模型可視化，這裏面咱們僅僅可視化基礎的雙高斯獨立分佈，代碼以下，

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
from matplotlib import cm
import matplotlib as mpl

num = 200
l = np.linspace(-5,5,num)
X, Y =np.meshgrid(l, l)
pos = np.concatenate((np.expand_dims(X,axis=2),np.expand_dims(Y,axis=2)),axis=2)

u1 = np.array([2, 2])
o1 = 3*np.array([[1, 0],
              [0, 1]])
a1 = (pos-u1).dot(np.linalg.inv(o1))
b1 = np.expand_dims(pos-u1,axis=3)
Z1 = np.zeros((num,num), dtype=np.float32)

u2 = np.array([-2, -2])
o2 = 3*np.array([[1, 0],
              [0, 1]])
a2 = (pos-u2).dot(np.linalg.inv(o2))
b2 = np.expand_dims(pos-u2,axis=3)
Z2 = np.zeros((num,num), dtype=np.float32)

for i in range(num):
    Z1[i] = [np.dot(a1[i,j],b1[i,j]) for j in range(num)]
    Z2[i] = [np.dot(a2[i,j],b2[i,j]) for j in range(num)]
Z1 = np.exp(Z1*(-0.5))/(2*np.pi*np.linalg.det(o1))
Z2 = np.exp(Z2*(-0.5))/(2*np.pi*np.linalg.det(o1))

Z = Z1 + Z2

fig = plt.figure()
ax = fig.add_subplot(211,projection='3d')
ax.plot_surface(X, Y, Z, rstride=5, cstride=5, alpha=0.3, cmap=cm.coolwarm)

cset = ax.contour(X,Y,Z,10,zdir='z',offset=0,cmap=cm.coolwarm)
cset = ax.contour(X, Y, Z, zdir='x', offset=-5,cmap=mpl.cm.winter)
cset = ax.contour(X, Y, Z, zdir='y', offset= 5,cmap= mpl.cm.winter)
'''
mpl.cm.rainbow
mpl.cm.winter
mpl.cm.bwr  # 藍，白，紅
cm.coolwarm
'''

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()

ax2 = fig.add_subplot(212)
cs = ax2.contour(X,Y,Z)
ax2.clabel(cs, inline=1, fontsize=20)

不過吳老大的圖兩個高斯分佈投影是分開的，因此咱們再次小改繪圖部分，

cset = ax.contour(X,Y,Z1,10,zdir='z',offset=0,cmap=cm.coolwarm)
cset = ax.contour(X,Y,Z2,10,zdir='z',offset=0,cmap=cm.coolwarm)
cset = ax.contour(X, Y, Z, zdir='x', offset=-5,cmap=mpl.cm.winter)
cset = ax.contour(X, Y, Z, zdir='y', offset= 5,cmap= mpl.cm.winter)
'''
mpl.cm.rainbow
mpl.cm.winter
mpl.cm.bwr  # 藍，白，紅
cm.coolwarm
'''

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()

ax2 = fig.add_subplot(212)
cs = ax2.contour(X,Y,Z1)
ax2.clabel(cs, inline=1, fontsize=20)
cs2 = ax2.contour(X,Y,Z2)
ax2.clabel(cs2, inline=1, fontsize=20)

顯示以下，框子不夠標準致使圓有點變形，不過這個能夠經過手動拉伸獲得優化，因此問題不大，

有關多元正態分佈的數學原理建議自行百度（cs229的學習不會在博客上更新，主要是由於我很是很是討厭打數學公式233）。