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創建

標準Python的列表(list)中，元素本質是對象。

如：L = [1, 2, 3]，需要3個指針和三個整數對象，對于數值運算比較浪費內存和CPU。

因此，Numpy提供了ndarray(N-dimensional array object)對象：存儲單一數據類型的多維數組。

# 通過array函數傳遞list對象

L = [1, 2, 3, 4, 5, 6]

print ("L = ", L)

a = np.array(L) # 數組沒有逗號

print ("a = ", a)

print (type(a), type(L))

# 若傳遞的是多層嵌套的list，將創建多維數組

b = np.array([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]])

print (b)

# 數組大小可以通過其shape屬性獲得

print (a.shape)

print (b.shape)

# # 也可以強制修改shape

b.shape = 4, 3

print (b)

# 注：從(3,4)改為(4,3)并不是對數組進行轉置，而只是改變每個軸的大小，數組元素在內存中的位置并沒有改變

# 當某個軸為-1時，將根據數組元素的個數自動計算此軸的長度

b.shape = 2, -1

print (b)

print (b.shape)

b.shape = 3, 4

print (b)

# 使用reshape方法，可以創建改變了尺寸的新數組，原數組的shape保持不變

c = b.reshape((4, -1))

print ("b = \n", b)

print ('c = \n', c)

# 數組b和c共享內存，修改任意一個將影響另外一個

b[0][1] = 20

print ("b = \n", b)

print ("c = \n", c)

# 數組的元素類型可以通過dtype屬性獲得

print (a.dtype)

print (b.dtype)

# 可以通過dtype參數在創建時指定元素類型

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安全的轉換

f = d.astype(np.int)

print (f)

# 但不要強制僅修改元素類型，如下面這句，將會以int來解釋單精度float類型

d.dtype = np.int

print (d)

2.使用函數創建

如果生成一定規則的數據，可以使用NumPy提供的專門函數

arange函數類似于python的range函數：指定起始值、終止值和步長來創建數組，和Python的range類似，arange同樣不包括終值；但arange可以生成浮點類型，而range只能是整數類型

a = np.arange(1, 10, 0.5)

print (a)

# linspace函數通過指定起始值、終止值和元素個數來創建數組，缺省包括終止值

b = np.linspace(1, 10, 10)

print ('b = ', b)

# 可以通過endpoint關鍵字指定是否包括終值

c = np.linspace(1, 10, 10, endpoint=False)

print ('c = ', c)

# 和linspace類似，logspace可以創建等比數列

# 下面函數創建起始值為10^1，終止值為10^2，有10個數的等比數列

d = np.logspace(1, 2, 9, endpoint=True)

print (d)

# 下面創建起始值為2^0，終止值為2^10(包括)，有10個數的等比數列

f = np.logspace(0, 10, 10, endpoint=True, base=2)

print (f)

# 使用 frombuffer, fromstring, fromfile等函數可以從字節序列創建數組

s = 'abcdz'

g = np.fromstring(s, dtype=np.int8) # 復制對應的ASII碼

print (g)

3.存取

3.1常規方法

# 數組元素的存取方法和Python的標準方法相同

a = np.arange(10)

print (a)

# 獲取某個元素

print (a[3])

# 切片[3,6)，左閉右開

print (a[3:6])

# 省略開始下標，表示從0開始

print (a[:5])

# 下標為負表示從后向前數

print (a[3:])

# 步長為2

print (a[1:9:2])

# 步長為-1，即翻轉

print (a[::-1])

# 切片數據是原數組的一個視圖，與原數組共享內容空間，可以直接修改元素值

a[1:4] = 10, 20, 30

print (a)

# 因此，在實踐中，切實注意原始數據是否被破壞，如：

b = a[2:5]

b[0] = 200

print (a)

3.2 整數/布爾數組存取

3.2.1

# 根據整數數組存取：當使用整數序列對數組元素進行存取時，

# 將使用整數序列中的每個元素作為下標，整數序列可以是列表(list)或者數組(ndarray)。

# 使用整數序列作為下標獲得的數組不和原始數組共享數據空間。

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

# 使用布爾數組i作為下標存取數組a中的元素：返回數組a中所有在數組b中對應下標為True的元素

# 生成10個滿足[0,1)中均勻分布的隨機數

a = np.random.rand(10)

print (a)

# 大于0.5的元素索引

print (a > 0.5)

# 大于0.5的元素

b = a[a > 0.5]

print (b)

# 將原數組中大于0.5的元素截取成0.5

a[a > 0.5] = 0.5

print (a)

# # # b不受影響

print (b)

3.3 二維數組的切片

a = np.arange(0, 60, 10) # 行向量

print ('a = ', a)

b = a.reshape((-1, 1)) # 轉換成列向量

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)

# 二維數組的切片

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數學庫的時間比較

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直接使用庫函數

a = np.array((1, 2, 3, 4, 5, 5, 7, 3, 2, 2, 8, 8))

print ('原始數組：', a)

# 使用庫函數unique

b = np.unique(a)

print ('去重后：', b)

4.2.2 二維數組的去重，結果會是預期的么？

c = np.array(((1, 2), (3, 4), (5, 6), (1, 3), (3, 4), (7, 6)))

print (u'二維數組：\n', c)

print ('去重后：', np.unique(c))

4.2.3 方案1：轉換為虛數

# r, i = np.split(c, (1, ), axis=1)

# x = r + i * 1j

x = c[:, 0] + c[:, 1] * 1j

print ('轉換成虛數：', x)

print ('虛數去重后：', np.unique(x))

print (np.unique(x, return_index=True)) # 思考return_index的意義

idx = np.unique(x, return_index=True)[1]

print ('二維數組去重：\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 繪制正態分布概率密度函數

# 自定義字體使中文正常顯示

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'高斯分布函數', fontsize=18)

plt.grid(True)

plt.show()

5.2 損失函數

# 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'])]

歡迎關注我的公眾號，新手入門資料。

wechat.jpg

總結

以上是生活随笔為你收集整理的自学python编程基础科学计算_Python基础与科学计算常用方法的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。