Python基礎與科學計算常用方法
本文使用的是Jupyter Notebook竹捉，Python3烫止。你可以將代碼直接復制到Jupyter Notebook中運行，以便更好的學習。
導入所需要的頭文件
import numpy as np
import numpy as np
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import time
from scipy.optimize import leastsq
from scipy import stats
import scipy.optimize as opt
import matplotlib.pyplot as plt
from scipy.stats import norm, poisson
import math
import scipy
from scipy.interpolate import BarycentricInterpolator
from scipy.interpolate import CubicSpline
a = np.arange(0, 60, 10).reshape((-1, 1)) + np.arange(6)
print (a)
1.使用array創(chuàng)建
標準Python的列表(list)中话侧，元素本質(zhì)是對象。
如：L = [1, 2, 3]闯参，需要3個指針和三個整數(shù)對象瞻鹏，對于數(shù)值運算比較浪費內(nèi)存和CPU。
因此鹿寨，Numpy提供了ndarray(N-dimensional array object)對象：存儲單一數(shù)據(jù)類型的多維數(shù)組新博。
# 通過array函數(shù)傳遞list對象
L = [1, 2, 3, 4, 5, 6]
print ("L = ", L)
a = np.array(L) # 數(shù)組沒有逗號
print ("a = ", a)
print (type(a), type(L))
# 若傳遞的是多層嵌套的list，將創(chuàng)建多維數(shù)組
b = np.array([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]])
print (b)
# 數(shù)組大小可以通過其shape屬性獲得
print (a.shape)
print (b.shape)
# # 也可以強制修改shape
b.shape = 4, 3
print (b)
# 注：從(3,4)改為(4,3)并不是對數(shù)組進行轉(zhuǎn)置脚草，而只是改變每個軸的大小赫悄，數(shù)組元素在內(nèi)存中的位置并沒有改變
# 當某個軸為-1時，將根據(jù)數(shù)組元素的個數(shù)自動計算此軸的長度
b.shape = 2, -1
print (b)
print (b.shape)
b.shape = 3, 4
print (b)
# 使用reshape方法馏慨，可以創(chuàng)建改變了尺寸的新數(shù)組埂淮，原數(shù)組的shape保持不變
c = b.reshape((4, -1))
print ("b = \n", b)
print ('c = \n', c)
# 數(shù)組b和c共享內(nèi)存，修改任意一個將影響另外一個
b[0][1] = 20
print ("b = \n", b)
print ("c = \n", c)
# 數(shù)組的元素類型可以通過dtype屬性獲得
print (a.dtype)
print (b.dtype)
# 可以通過dtype參數(shù)在創(chuàng)建時指定元素類型
d = np.array([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]], dtype=np.float)
f = np.array([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]], dtype=np.complex)
print (d)
print (f)
# 如果更改元素類型写隶，可以使用astype安全的轉(zhuǎn)換
f = d.astype(np.int)
print (f)
# 但不要強制僅修改元素類型倔撞，如下面這句，將會以int來解釋單精度float類型
d.dtype = np.int
print (d)
2.使用函數(shù)創(chuàng)建
如果生成一定規(guī)則的數(shù)據(jù)慕趴，可以使用NumPy提供的專門函數(shù)
arange函數(shù)類似于python的range函數(shù)：指定起始值痪蝇、終止值和步長來創(chuàng)建數(shù)組，和Python的range類似冕房，arange同樣不包括終值躏啰；但arange可以生成浮點類型，而range只能是整數(shù)類型
a = np.arange(1, 10, 0.5)
print (a)
# linspace函數(shù)通過指定起始值耙册、終止值和元素個數(shù)來創(chuàng)建數(shù)組丙唧，缺省包括終止值
b = np.linspace(1, 10, 10)
print ('b = ', b)
# 可以通過endpoint關(guān)鍵字指定是否包括終值
c = np.linspace(1, 10, 10, endpoint=False)
print ('c = ', c)
# 和linspace類似，logspace可以創(chuàng)建等比數(shù)列
# 下面函數(shù)創(chuàng)建起始值為10^1觅玻，終止值為10^2想际，有10個數(shù)的等比數(shù)列
d = np.logspace(1, 2, 9, endpoint=True)
print (d)
# 下面創(chuàng)建起始值為2^0，終止值為2^10(包括)溪厘，有10個數(shù)的等比數(shù)列
f = np.logspace(0, 10, 10, endpoint=True, base=2)
print (f)
# 使用 frombuffer, fromstring, fromfile等函數(shù)可以從字節(jié)序列創(chuàng)建數(shù)組
s = 'abcdz'
g = np.fromstring(s, dtype=np.int8) # 復制對應的ASII碼
print (g)
3.存取
3.1常規(guī)方法
# 數(shù)組元素的存取方法和Python的標準方法相同
a = np.arange(10)
print (a)
# 獲取某個元素
print (a[3])
# 切片[3,6)胡本，左閉右開
print (a[3:6])
# 省略開始下標，表示從0開始
print (a[:5])
# 下標為負表示從后向前數(shù)
print (a[3:])
# 步長為2
print (a[1:9:2])
# 步長為-1畸悬，即翻轉(zhuǎn)
print (a[::-1])
# 切片數(shù)據(jù)是原數(shù)組的一個視圖侧甫，與原數(shù)組共享內(nèi)容空間，可以直接修改元素值
a[1:4] = 10, 20, 30
print (a)
# 因此，在實踐中披粟，切實注意原始數(shù)據(jù)是否被破壞咒锻，如：
b = a[2:5]
b[0] = 200
print (a)
3.2 整數(shù)/布爾數(shù)組存取
3.2.1
# 根據(jù)整數(shù)數(shù)組存取：當使用整數(shù)序列對數(shù)組元素進行存取時守屉，
# 將使用整數(shù)序列中的每個元素作為下標惑艇，整數(shù)序列可以是列表(list)或者數(shù)組(ndarray)。
# 使用整數(shù)序列作為下標獲得的數(shù)組不和原始數(shù)組共享數(shù)據(jù)空間拇泛。
a = np.logspace(0, 9, 10, base=2)
print (a)
i = np.arange(0, 10, 2)
print (i)
# 利用i取a中的元素
b = a[i]
print (b)
# b的元素更改滨巴，a中元素不受影響
b[2] = 1.6
print (b)
print (a)
3.2.2
# 使用布爾數(shù)組i作為下標存取數(shù)組a中的元素：返回數(shù)組a中所有在數(shù)組b中對應下標為True的元素
# 生成10個滿足[0,1)中均勻分布的隨機數(shù)
a = np.random.rand(10)
print (a)
# 大于0.5的元素索引
print (a > 0.5)
# 大于0.5的元素
b = a[a > 0.5]
print (b)
# 將原數(shù)組中大于0.5的元素截取成0.5
a[a > 0.5] = 0.5
print (a)
# # # b不受影響
print (b)
3.3 二維數(shù)組的切片
a = np.arange(0, 60, 10) # 行向量
print ('a = ', a)
b = a.reshape((-1, 1)) # 轉(zhuǎn)換成列向量
print (b)
c = np.arange(6)
print (c)
f = b + c # 行 + 列
print (f)
# 合并上述代碼：
a = np.arange(0, 60, 10).reshape((-1, 1)) + np.arange(6)
print (a)
# 二維數(shù)組的切片
print (a[[0, 1, 2], [2, 3, 4]])
print (a[4, [2, 3, 4]])
print (a[4:, [2, 3, 4]])
i = np.array([True, False, True, False, False, True])
print (a[i])
print (a[i, 3])
4.1 numpy與Python數(shù)學庫的時間比較
for j in np.logspace(0, 7, 10):
j = int(j)
x = np.linspace(0, 10, j)
start = time.clock()
y = np.sin(x)
t1 = time.clock() - start
x = x.tolist()
start = time.clock()
for i, t in enumerate(x):
x[i] = math.sin(t)
t2 = time.clock() - start
print (j, ": ", t1, t2, t2/t1)
# 4.2 元素去重
4.2.1直接使用庫函數(shù)
a = np.array((1, 2, 3, 4, 5, 5, 7, 3, 2, 2, 8, 8))
print ('原始數(shù)組：', a)
# 使用庫函數(shù)unique
b = np.unique(a)
print ('去重后：', b)
4.2.2 二維數(shù)組的去重，結(jié)果會是預期的么俺叭？
c = np.array(((1, 2), (3, 4), (5, 6), (1, 3), (3, 4), (7, 6)))
print (u'二維數(shù)組：\n', c)
print ('去重后：', np.unique(c))
4.2.3 方案1：轉(zhuǎn)換為虛數(shù)
# r, i = np.split(c, (1, ), axis=1)
# x = r + i * 1j
x = c[:, 0] + c[:, 1] * 1j
print ('轉(zhuǎn)換成虛數(shù)：', x)
print ('虛數(shù)去重后：', np.unique(x))
print (np.unique(x, return_index=True)) # 思考return_index的意義
idx = np.unique(x, return_index=True)[1]
print ('二維數(shù)組去重：\n', c[idx])
4.2.3 方案2：利用set
print ('去重方案2：\n', np.array(list(set([tuple(t) for t in c]))))
4.3 stack and axis
a = np.arange(1, 10).reshape((3, 3))
b = np.arange(11, 20).reshape((3, 3))
c = np.arange(101, 110).reshape((3, 3))
print ('a = \n', a)
print ('b = \n', b)
print ('c = \n', c)
print ('axis = 0 \n', np.stack((a, b, c), axis=0))
print ('axis = 1 \n', np.stack((a, b, c), axis=1))
print ('axis = 2 \n', np.stack((a, b, c), axis=2))
a = np.arange(1, 10).reshape(3,3)
print (a)
b = a + 10
print (b)
print (np.dot(a, b)) # dot 正常的矩陣乘法
print (a * b) # * 對應元素的相乘
a = np.arange(1, 10)
print (a)
b = np.arange(20,25)
print (b)
print (np.concatenate((a, b)))
5.繪圖
5.1 繪制正態(tài)分布概率密度函數(shù)
# 自定義字體使中文正常顯示
mpl.rcParams['font.sans-serif'] = [u'SimHei'] #FangSong/黑體 FangSong/KaiTi
mpl.rcParams['axes.unicode_minus'] = False
mu = 0
sigma = 1
x = np.linspace(mu - 3 * sigma, mu + 3 * sigma, 51)
y = np.exp(-(x - mu) ** 2 / (2 * sigma ** 2)) / (math.sqrt(2 * math.pi) * sigma)
print (x.shape)
print ('x = \n', x)
print (y.shape)
print ('y = \n', y)
plt.figure(facecolor='w') # 背景色設置為白色
# plt.plot(x, y, 'ro-', linewidth=2)
plt.plot(x, y, 'r-', x, y, 'go', linewidth=2, markersize=8)
plt.xlabel('X', fontsize=15)
plt.ylabel('Y', fontsize=15)
plt.title(u'高斯分布函數(shù)', fontsize=18)
plt.grid(True)
plt.show()
5.2 損失函數(shù)
# Logistic損失(-1,1)/SVM Hinge損失/ 0/1損失
plt.figure(figsize=(10,8),dpi=100) # 指定圖像尺寸和dpi
x = np.array(np.linspace(start=-2, stop=3, num=1001, dtype=np.float))
y_logit = np.log(1 + np.exp(-x)) / math.log(2)
y_boost = np.exp(-x)
y_01 = x < 0
y_hinge = 1.0 - x
y_hinge[y_hinge < 0] = 0
plt.plot(x, y_logit, 'r-', label='Logistic Loss', linewidth=2)
plt.plot(x, y_01, 'g-', label='0/1 Loss', linewidth=2)
plt.plot(x, y_hinge, 'b-', label='Hinge Loss', linewidth=2)
plt.plot(x, y_boost, 'm--', label='Adaboost Loss', linewidth=2)
plt.grid() # 畫格子出來
plt.legend(loc='upper right') # 圖例顯示位置
plt.savefig('1.png')
plt.show()
5.3 x^x
def f(x):
y = np.ones_like(x)
i = x > 0
y[i] = np.power(x[i], x[i])
i = x < 0
y[i] = np.power(-x[i], -x[i])
return y
x = np.linspace(-1.3, 1.3, 101)
y = f(x)
plt.plot(x, y, 'g-', label='x^x', linewidth=2)
plt.grid()
plt.legend(loc='upper left')
plt.show()
5.4 胸型線
x = np.arange(1, 0, -0.001)
y = (-3 * x * np.log(x) + np.exp(-(40 * (x - 1 / np.e)) ** 4) / 25) / 2
plt.figure(figsize=(5,7), facecolor='w')
plt.plot(y, x, 'r-', linewidth=2)
plt.grid(True)
plt.title(u'胸型線', fontsize=20)
# plt.savefig('breast.png')
plt.show()
5.5 心形線
t = np.linspace(0, 2*np.pi, 100)
x = 16 * np.sin(t) ** 3
y = 13 * np.cos(t) - 5 * np.cos(2*t) - 2 * np.cos(3*t) - np.cos(4*t)
plt.plot(x, y, 'r-', linewidth=2)
plt.grid(True)
plt.show()
5.6 漸開線
t = np.linspace(0, 50, num=1000)
x = t*np.sin(t) + np.cos(t)
y = np.sin(t) - t*np.cos(t)
plt.plot(x, y, 'r-', linewidth=2)
plt.grid()
plt.show()
5.7Bar
x = np.arange(0, 10, 0.1)
y = np.sin(x)
plt.bar(x, y, width=0.04, linewidth=0.2)
plt.plot(x, y, 'r--', linewidth=2)
plt.title(u'Sin曲線')
plt.xticks(rotation=-60)
plt.xlabel('X')
plt.ylabel('Y')
plt.grid()
plt.show()
6. 概率分布
6.1 均勻分布
x = np.random.rand(10000)
t = np.arange(len(x))
# plt.hist(x, 30, color='#000000', alpha=0.5, label=u'均勻分布') # alpha透明度
plt.plot(t, x, 'g.', label=u'均勻分布') # .小點, o 圈, - 線
plt.legend(loc='upper left')
plt.grid()
plt.show()
6.2 驗證中心極限定理
t = 1000
a = np.zeros(10000)
for i in range(t):
a += np.random.uniform(-5, 5, 10000)
a /= t
plt.hist(a, bins=30, color='g', alpha=0.5, normed=True, label=u'均勻分布疊加')
plt.legend(loc='upper left')
plt.grid()
plt.show()
# 6.21 其他分布的中心極限定理
lamda = 10
p = stats.poisson(lamda)
y = p.rvs(size=1000)
mx = 30
r = (0, mx)
bins = r[1] - r[0]
plt.figure(figsize=(10, 8), facecolor='w')
plt.subplot(121)
plt.hist(y, bins=bins, range=r, color='g', alpha=0.8, normed=True)
t = np.arange(0, mx+1)
plt.plot(t, p.pmf(t), 'ro-', lw=2)
plt.grid(True)
N = 1000
M = 10000
plt.subplot(122)
a = np.zeros(M, dtype=np.float)
p = stats.poisson(lamda)
for i in np.arange(N):
y = p.rvs(size=M)
a += y
a /= N
plt.hist(a, bins=20, color='g', alpha=0.8, normed=True)
plt.grid(b=True)
plt.show()
6.3 Poisson分布
x = np.random.poisson(lam=5, size=10000)
print (x)
pillar = 15
a = plt.hist(x, bins=pillar, normed=True, range=[0, pillar], color='g', alpha=0.5)
plt.grid()
# plt.show()
print (a)
print (a[0].sum())
6.4 直方圖的使用
mu = 2
sigma = 3
data = mu + sigma * np.random.randn(1000)
h = plt.hist(data, 30, normed=1, color='#a0a0ff')
x = h[1]
y = norm.pdf(x, loc=mu, scale=sigma)
plt.plot(x, y, 'r--', x, y, 'ro', linewidth=2, markersize=4)
plt.grid()
plt.show()
6.5 插值
rv = poisson(5)
x1 = a[1]
y1 = rv.pmf(x1)
itp = BarycentricInterpolator(x1, y1) # 重心插值
x2 = np.linspace(x.min(), x.max(), 50)
y2 = itp(x2)
cs = scipy.interpolate.CubicSpline(x1, y1) # 三次樣條插值
plt.plot(x2, cs(x2), 'm--', linewidth=5, label='CubicSpine') # 三次樣條插值
plt.plot(x2, y2, 'g-', linewidth=3, label='BarycentricInterpolator') # 重心插值
plt.plot(x1, y1, 'r-', linewidth=1, label='Actural Value') # 原始值
plt.legend(loc='upper right')
plt.grid()
plt.show()
7. 繪制三維圖像
# x, y = np.ogrid[-3:3:100j, -3:3:100j]
# print(x,y)
u = np.linspace(-3, 3, 101)
x, y = np.meshgrid(u, u)
z = x*y*np.exp(-(x**2 + y**2)/2) / math.sqrt(2*math.pi)
# z = x*y*np.exp(-(x**2 + y**2)/2) / math.sqrt(2*math.pi)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
# ax.plot_surface(x, y, z, rstride=5, cstride=5, cmap=cm.coolwarm, linewidth=0.1) #
ax.plot_surface(x, y, z, rstride=5, cstride=5, cmap=cm.Accent, linewidth=0.5)
plt.show()
# cmaps = [('Perceptually Uniform Sequential',
# ['viridis', 'inferno', 'plasma', 'magma']),
# ('Sequential', ['Blues', 'BuGn', 'BuPu',
# 'GnBu', 'Greens', 'Greys', 'Oranges', 'OrRd',
# 'PuBu', 'PuBuGn', 'PuRd', 'Purples', 'RdPu',
# 'Reds', 'YlGn', 'YlGnBu', 'YlOrBr', 'YlOrRd']),
# ('Sequential (2)', ['afmhot', 'autumn', 'bone', 'cool',
# 'copper', 'gist_heat', 'gray', 'hot',
# 'pink', 'spring', 'summer', 'winter']),
# ('Diverging', ['BrBG', 'bwr', 'coolwarm', 'PiYG', 'PRGn', 'PuOr',
# 'RdBu', 'RdGy', 'RdYlBu', 'RdYlGn', 'Spectral',
# 'seismic']),
# ('Qualitative', ['Accent', 'Dark2', 'Paired', 'Pastel1',
# 'Pastel2', 'Set1', 'Set2', 'Set3']),
# ('Miscellaneous', ['gist_earth', 'terrain', 'ocean', 'gist_stern',
# 'brg', 'CMRmap', 'cubehelix',
# 'gnuplot', 'gnuplot2', 'gist_ncar',
# 'nipy_spectral', 'jet', 'rainbow',
# 'gist_rainbow', 'hsv', 'flag', 'prism'])]
歡迎關(guān)注我的公眾號恭取，閱讀更多新手入門資料。
wechat.jpg
