矩陣
矩阵是numpy.matrix类类型的对象，该类继承自numpy.ndarray，任何针对多维数组的操作，对矩阵同样有效，但是作为子类矩阵又结合其自身的特点，做了必要的扩充，比如：乘法计算、求逆等。
- 矩阵对象的创建可以通过以下三种方式：
numpy.matrix(任何可被解释为矩阵的二维容器,copy=是否复制数据(缺省值为True，即复制数据))- 如果copy的值为True(缺省)，所得到的矩阵对象与参数中的源容器共享同一份数据，否则，各自拥有独立的数据拷贝。
numpy.mat(任何可被解释为矩阵的二维容器)等价于numpy.matrix(..., copy=False)- 由该函数创建的矩阵对象与参数中的源容器一定共享数据，无法拥有独立的数据拷贝。
numpy.bmat(拼块规则)- 包含若干小矩阵数据的大矩阵，拼接规则由参数指定
以上函数也可以接受字符串形式的矩阵描述：数据项通过空格分隔，数据行通过分号分隔。例如：
'1 2 3; 4 5 6' / 1 2 3 4 5 6 /
#mat.pyimport numpy as npa = np.array([ [1, 2, 6], [3, 5, 7], [4, 8, 9]])print(a, type(a))b = np.matrix(a)print(b, type(b))b += 10print(b)print(a)c = np.matrix(a, copy=False)print(c, type(c))c += 10print(c)print(a)d = np.mat(a)print(d, type(d))d -= 10print(d)print(a)e = np.mat('1 2 6; 3 5 7; 4 8 9')print(e)f = np.bmat('b e')print(f)g = np.bmat('b e; e b')# 多維數組的乘法：對應元素相乘h = a * aprint('h等於\n', h)# 矩陣相乘：乘積矩陣的第i行第j列元素等於被乘矩陣的第i行於乘數矩陣的第j列的點積i = e * eprint('i等於\n', i)j = e.Iprint('e的逆矩陣為\n', j)print(j * e)k = a.dot(a) # 數組使用矩陣乘法l = np.linalg.inv(a) # 數組使用逆矩陣 [[1 2 6] [3 5 7] [4 8 9]][[1 2 6] [3 5 7] [4 8 9]] [[11 12 16] [13 15 17] [14 18 19]][[1 2 6] [3 5 7] [4 8 9]][[1 2 6] [3 5 7] [4 8 9]] [[11 12 16] [13 15 17] [14 18 19]][[11 12 16] [13 15 17] [14 18 19]][[11 12 16] [13 15 17] [14 18 19]] [[1 2 6] [3 5 7] [4 8 9]][[1 2 6] [3 5 7] [4 8 9]][[1 2 6] [3 5 7] [4 8 9]][[11 12 16 1 2 6] [13 15 17 3 5 7] [14 18 19 4 8 9]]h等於 [[ 1 4 36] [ 9 25 49] [16 64 81]]i等於 [[ 31 60 74] [ 46 87 116] [ 64 120 161]]e的逆矩陣為 [[-0.73333333 2. -1.06666667] [ 0.06666667 -1. 0.73333333] [ 0.26666667 0. -0.06666667]][[ 1.00000000e+00 1.77635684e-15 3.55271368e-15] [ 0.00000000e+00 1.00000000e+00 -1.11022302e-16] [ 0.00000000e+00 0.00000000e+00 1.00000000e+00]]
通用函数
1. frompyfunc->ufunc对象
def 标量函数(标量参数1, 标量参数2, ...): ... return 标量返回值1, 标量返回值2, ...矢量参数1矢量参数2...numpy.frompyfunc(标量函数, 参数个数, 返回值个数) ->矢量函数 # numpy.ufunc类类型的对象，可调用对象矢量函数(矢量参数1, 矢量参数2, ...) ->矢量返回值1, 矢量返回值2
# ufunc.pyimport numpy as npdef foo(x, y): return x + y, x - y, x * ydef hum(x): def fun(y): return x + y, x - y, x * y return np.frompyfunc(fun, 1, 1)x, y = 1, 4print(foo(x, y))X, Y = np.array([1, 2, 3]), np.array([4, 5, 6])bar = np.frompyfunc(foo, 2, 3)print(bar(X, Y))print(np.frompyfunc(foo, 2, 3)(X, Y))print(hum(100)(X))
(5, -3, 4)(array([5, 7, 9], dtype=object), array([-3, -3, -3], dtype=object), array([4, 10, 18], dtype=object))(array([5, 7, 9], dtype=object), array([-3, -3, -3], dtype=object), array([4, 10, 18], dtype=object))[(101, 99, 100) (102, 98, 200) (103, 97, 300)]
2. 加法通用函数：add
add.reduce() - 累加和add.accumulate() - 累加和过程add.reduceat() - 分段累加和add.outer() - 外和
# add.pya = np.arange(1, 7)print(a)b = a + aprint(b)b = np.add(a, a)print(b)c = np.add.reduce(a)print(c)d = np.add.accumulate(a)print(d)e = np.add.reduceat(a, [0, 2, 4]) # 按照下標分段print(e)f = np.add.outer([10, 20, 30], a)print(f)g = np.outer([10, 20, 30], a) # 不加add的話是外積print(g)
[1 2 3 4 5 6][ 2 4 6 8 10 12][ 2 4 6 8 10 12]21[ 1 3 6 10 15 21][ 3 7 11][[11 12 13 14 15 16] [21 22 23 24 25 26] [31 32 33 34 35 36]][[ 10 20 30 40 50 60] [ 20 40 60 80 100 120] [ 30 60 90 120 150 180]]
3. 除法通用函数
[5 5 -5 -5] <真除> [2 -2 2 -2] = [2.5 -2.5 -2.5 2.5] numpy.true_divide() numpy.divide() / [5 5 -5 -5] <地板除> [2 -2 2 -2] = [2 -3 -3 2] numpy.floor_divide() // [5 5 -5 -5] <天花板除> [2 -2 2 -2] = [3 -2 -2 3] 天花板取整(真除的结果)：numpy.ceil() [5 5 -5 -5] <截断除> [2 -2 2 -2] = [2 -2 -2 2] 截断取整(真除的结果)：numpy.trunc()
# div.pya = np.array([5, 5, -5, -5])b = np.array([2, -2, 2, -2])print(a, b)# c = np.true_divide(a, b)# c = np.divide(a, b)c = a / bprint(c)# d = np.floor_divide(a, b)d = a // bprint(d)e = np.ceil(a / b).astype(int)print(e)f = np.trunc(a / b).astype(int)print(f)
[ 5 5 -5 -5] [ 2 -2 2 -2][ 2.5 -2.5 -2.5 2.5][ 2 -3 -3 2][ 3 -2 -2 3][ 2 -2 -2 2]
4. 取余通用函数
被除数/除数=商...余数除数 * 商 + 余数 = 被除数[5 5 -5 -5] <地板除> [2 -2 2 -2] = [2 -3 -3 2]...[1 -1 1 -1]numpy.remainder()numpy.mod()%[5 5 -5 -5] <截断除> [2 -2 2 -2] = [2 -2 -2 2]...[1 1 -1 -1]numpy.fmod()
Numpy将Python语言中针对标量的运算符，通过通用函数加以重载定义，以支持数组形式的矢量运算。
# fib.pyimport numpy as npn = 35# 递归方法# def fibo(n):# return 1 if n < 3 else fibo(n - 1) + fibo(n - 2)# print(fibo(n))# 循环换位方法fn_1, fn_2 = 0, 1for i in range(n): fn = fn_1 + fn_2 fn_1, fn_2 = fn, fn_1print(fn)# 矩阵乘方print(int((np.mat('1. 1.; 1. 0.') ** (n - 1))[0, 0]))# 斐波那契公式r = np.sqrt(5)print(int((((1 + r) / 2) ** n - ((1 - r) / 2) ** n) / r)) 922746592274659227465
5. Numpy中的所有三角函数都是通用函数
x = Asin(at+pi/2)y = Bsin(bt)
# lissa.pyimport numpy as npimport matplotlib.pyplot as mpt = np.linspace(0, 2 * np.pi, 201)A, a, B, b = 10, 9, 5, 8x = A * np.sin(a * t + np.pi / 2)y = B * np.sin(b * t)mp.figure('Lissajous', facecolor='lightgray')mp.title('Lissajous', fontsize=20)mp.xlabel('x', fontsize=14)mp.ylabel('y', fontsize=14)mp.tick_params(labelsize=10)mp.grid(linestyle=':')mp.plot(x, y, c='orangered', label='Lissajous')mp.legend()mp.show() 
y = 4/(1pi) sin(1x) 2 x 1 - 1y = 4/(3pi) sin(3x) 2 x 2 - 1y = 4/(5pi) sin(5x) 2 x 3 - 1
import numpy as npimport matplotlib.pyplot as mpdef squarewave(n): k = np.arange(1, n + 1) def fun(x): return np.sum(4 / ((2 * k - 1) * np.pi) * np.sin((2 * k - 1) * x)) return np.frompyfunc(fun, 1, 1)x = np.linspace(0, 2 * np.pi, 201)y1 = squarewave(1)(x)y2 = squarewave(2)(x)y3 = squarewave(3)(x)y4 = squarewave(10)(x)y5 = squarewave(100)(x)y6 = squarewave(1000)(x)mp.figure('Squarewave', facecolor='lightgray')mp.title('Squarewave', fontsize=20)mp.xlabel('x', fontsize=14)mp.ylabel('y', fontsize=14)mp.tick_params(labelsize=10)mp.grid(linestyle=':')mp.plot(x, y1, label='n=1')mp.plot(x, y2, label='n=2')mp.plot(x, y3, label='n=3')mp.plot(x, y4, label='n=10')mp.plot(x, y5, label='n=100')mp.plot(x, y6, label='n=1000')mp.legend()mp.show() 
7. 位运算通用函数
位异或：^/__xor__/bitwise_xor0 ^ 0 = 00 ^ 1 = 11 ^ 0 = 11 ^ 1 = 0if a^b<0 then a和b异号位与：&/__and__/bitwise_and0 & 0 = 00 & 1 = 01 & 0 = 01 & 1 = 1# 判断是否为2的幂1 2^0 00001 0 000002 2^1 00010 1 000014 2^2 00100 3 000118 2^3 01000 7 0011116 2^4 10000 15 01111...if a & (a-1) == 0 then a是2的幂位或：|/__or__/bitwise_or0 | 0 = 00 | 1 = 11 | 0 = 11 | 1 = 1位反：~/__not__/bitwise_not~0 = 1~1 = 0移位：< >/__rshift__/right_shift左移1位相当于乘2，右移1位相当于除2。
import numpy as npa = np.array([0, -1, 2, -3, 4, -5])b = np.array([0, 1, 2, 3, 4, 5])print(a, b)# c = a ^ b# c = a.__xor__(b)c = np.bitwise_xor(a, b)print(c)print(np.where(c < 0)[0])d = np.arange(1, 21)print(d)# e = d & (d - 1)# e = d.__and__(d - 1)e = np.bitwise_and(d, d - 1)print(e)print(d[e == 0])# f = d << 1# f = d.__lshift__(1)f = np.left_shift(d, 1)print(f)# g = d >> 1# g = d.__rshift__(1)g = np.right_shift(d, 1)print(g)
[ 0 -1 2 -3 4 -5] [0 1 2 3 4 5][ 0 -2 0 -2 0 -2][1 3 5][ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20][ 0 0 2 0 4 4 6 0 8 8 10 8 12 12 14 0 16 16 18 16][ 1 2 4 8 16][ 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40][ 0 1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 9 9 10]
线性代数模块(linalg)
1. 逆矩阵和广义逆矩阵
如果一个方阵A与另一个方阵B的乘积是一个单位矩阵，那么A和B就互为逆矩阵。
np.linalg.inv(A)->A^-1将以上有关逆矩阵的定义推广到非方阵，则称为广义逆矩阵。
np.linalg.pinv(A)->A^-1 np.matrix.I -> inv/pinv
import numpy as npA = np.mat('1 2 3; 8 9 4; 7 6 5')print(A)B = np.linalg.inv(A)print(B)print(A * B)print(A.I)C = np.mat('11 12 13 14; 20 21 22 15; 19 18 17 16')print(C)# D = np.linalg.inv(C)D = np.linalg.pinv(C)print(D)print(C * D)print(C.I)E = np.mat('1 2 3; 4 5 6; 7 8 9')print(E)F = np.linalg.inv(E)print(F)print(E * F) [[1 2 3] [8 9 4] [7 6 5]][[-0.4375 -0.16666667 0.39583333] [ 0.25 0.33333333 -0.41666667] [ 0.3125 -0.16666667 0.14583333]][[ 1.00000000e+00 2.77555756e-17 -5.55111512e-17] [ 0.00000000e+00 1.00000000e+00 2.22044605e-16] [ 0.00000000e+00 1.94289029e-16 1.00000000e+00]][[-0.4375 -0.16666667 0.39583333] [ 0.25 0.33333333 -0.41666667] [ 0.3125 -0.16666667 0.14583333]][[11 12 13 14] [20 21 22 15] [19 18 17 16]][[-0.18055556 -0.08333333 0.23611111] [-0.04305556 0.04166667 -0.00138889] [ 0.09444444 0.16666667 -0.23888889] [ 0.1625 -0.125 0.0375 ]][[ 1.00000000e+00 2.22044605e-16 -3.33066907e-16] [ 1.77635684e-15 1.00000000e+00 -1.88737914e-15] [ 1.77635684e-15 8.88178420e-16 1.00000000e+00]][[-0.18055556 -0.08333333 0.23611111] [-0.04305556 0.04166667 -0.00138889] [ 0.09444444 0.16666667 -0.23888889] [ 0.1625 -0.125 0.0375 ]][[1 2 3] [4 5 6] [7 8 9]][[ 3.15251974e+15 -6.30503948e+15 3.15251974e+15] [-6.30503948e+15 1.26100790e+16 -6.30503948e+15] [ 3.15251974e+15 -6.30503948e+15 3.15251974e+15]][[ 0. 1. -0.5] [ 0. 2. -1. ] [ 0. 3. 2.5]]
2. 解线性方程组
$$
\left( \begin{array}{l} x - 2y + z = 0\ 2y - 8z - 8 = 0\- 4x + 5y + 9z + 9 = 0 \end{array} \right. $$
$$
\left( \begin{array}{l} 1x + - 2y + 1z = 0\ 0x + 2y + - 8z = 8\- 4x + 5y + 9z = - 9 \end{array} \right. $$
$$
\left( {\begin{array}{ {20}{c}} 1&{ - 2}&1\ 0&2&{ - 8}\ { - 4}&5&9 \end{array}} \right) \times \left( {\begin{array}{ {20}{c}} x\ y\ z \end{array}} \right) = \left( {\begin{array}{*{20}{c}} 0\ 8\ { - 9} \end{array}} \right) $$a x b = np.linalg.lstsq(a, b)[0] -> 拟合 = np.linalg.solve(a, b) -> 解
import numpy as npa = np.mat('1 -2 1; 0 2 -8; -4 5 9')b = np.mat('0; 8; -9')x = np.linalg.solve(a, b)print(x)x = np.linalg.lstsq(a, b)[0]print(x) [[29.] [16.] [ 3.]][[29.] [16.] [ 3.]]/Users/haoen110/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:6: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.
3. 特征值和特征向量
- 对于n阶方阵A，如果存在数a和非零n维列向量x，使得Ax=ax，则称a是矩阵A的一个特征值，x是矩阵A属于特征值a的特征向量 np.linalg.eig(A) -> 特征值数组，特征向量数组 [a1 a2] [[x11 x12] [x21 x22]]
import numpy as npA = np.mat('3 -2; 1 0')print(A)eigvals, eigvecs = np.linalg.eig(A)print(eigvals)print(eigvecs)print(A * eigvecs[:, 0])print(eigvals[0] * eigvecs[:, 0])print(A * eigvecs[:, 1])print(eigvals[1] * eigvecs[:, 1]) [[ 3 -2] [ 1 0]][2. 1.][[0.89442719 0.70710678] [0.4472136 0.70710678]][[1.78885438] [0.89442719]][[1.78885438] [0.89442719]][[0.70710678] [0.70710678]][[0.70710678] [0.70710678]]
4. 奇异值分解
- 主对角线上的元素称为矩阵M的奇异值，其它元素均为0 | M = U x S x V | | 正交矩阵 UxU^T = E = VxV^T np.linalg.svd(M, full_matrices=False)-> U, 奇异值, V
import numpy as npM = np.mat('4 11 14; 8 7 -2')print(M)U, sv, V = np.linalg.svd(M, full_matrices=False)print(U * U.T)print(V * V.T)print(sv)S = np.diag(sv)print(S)print(U * S * V) [[ 4 11 14] [ 8 7 -2]][[1.0000000e+00 3.2123061e-17] [3.2123061e-17 1.0000000e+00]][[ 1.00000000e+00 -6.16790569e-18] [-6.16790569e-18 1.00000000e+00]][18.97366596 9.48683298][[18.97366596 0. ] [ 0. 9.48683298]][[ 4. 11. 14.] [ 8. 7. -2.]]
5. 行列式
np.det(方阵)->行列式的值
import numpy as npA = np.mat('2 1; 3 4')print(A)print(np.linalg.det(A))B = np.mat('3 2 1; 4 9 8; 5 6 7')print(B)print(np.linalg.det(B)) [[2 1] [3 4]]5.000000000000001[[3 2 1] [4 9 8] [5 6 7]]47.999999999999986
快速傅里叶变换模块(fft)
- 原函数：y = f(x) - 时间/空间域函数
一系列正弦函数的叠加
y = A1sin(w1x+fai1) + A2sin(w2x+fai2) + ... + Ansin(wnx+fain) + R n->oo: R->0 [x1, x2, ..., xn]->[y1, y2, ..., yn]w1->A1, fai1 \ w2->A2, fai2 | A,fai = f'(w) - 频率域函数 ... | wn->An, fain / f(x) -傅里叶变换-> f'(w) 时空域 频率域 f(x) <-反傅里叶变换- f'(w) 时空域 频率域 np.fft.fftfreq(采样数, 采样周期)->频率序列 np.fft.fft(原函数值序列) -> 目标函数值序列(复数) 复数的模反映了振幅A，辐角反映了初相位fai np.fft.ifft(目标函数值序列(复数))->原函数值序列
import numpy as npimport numpy.fft as nfimport matplotlib.pyplot as mptimes = np.linspace(0, 2 * np.pi, 201)sigs1 = 4 / (1 * np.pi) * np.sin(1 * times)sigs2 = 4 / (3 * np.pi) * np.sin(3 * times)sigs3 = 4 / (5 * np.pi) * np.sin(5 * times)sigs4 = 4 / (7 * np.pi) * np.sin(7 * times)sigs5 = 4 / (9 * np.pi) * np.sin(9 * times)sigs6 = sigs1 + sigs2 + sigs3 + sigs4 + sigs5freqs = nf.fftfreq(times.size, times[1] - times[0])ffts = nf.fft(sigs6)pows = np.abs(ffts)sigs7 = nf.ifft(ffts).realmp.subplot(121)mp.title('Time Domain', fontsize=16)mp.xlabel('Time', fontsize=12)mp.ylabel('Signal', fontsize=12)mp.tick_params(labelsize=10)mp.grid(linestyle=':')mp.plot(times, sigs1, label='{:.4f}'.format( 1 / (2 * np.pi)))mp.plot(times, sigs2, label='{:.4f}'.format( 3 / (2 * np.pi)))mp.plot(times, sigs3, label='{:.4f}'.format( 5 / (2 * np.pi)))mp.plot(times, sigs4, label='{:.4f}'.format( 7 / (2 * np.pi)))mp.plot(times, sigs5, label='{:.4f}'.format( 9 / (2 * np.pi)))mp.plot(times, sigs6, label='{:.4f}'.format( 1 / (2 * np.pi)))mp.plot(times, sigs7, label='{:.4f}'.format( 1 / (2 * np.pi)), alpha=0.5, linewidth=6)mp.legend()mp.subplot(122)mp.title('Frequency Domain', fontsize=16)mp.xlabel('Frequency', fontsize=12)mp.ylabel('Power', fontsize=12)mp.tick_params(labelsize=10)mp.grid(linestyle=':')mp.plot(freqs[freqs >= 0], pows[freqs >= 0], c='orangered', label='Frequency Spectrum')mp.legend()mp.tight_layout() 
- 基于傅里叶变换的频域滤波 ____________________IFFT_____________________ | | v | 高能信号 FFT 频域滤波 | >含噪信号----->含噪频谱---------->高能频谱 低能噪声/
import numpy as npimport numpy.fft as nfimport scipy.io.wavfile as wfimport matplotlib.pyplot as mpsample_rate, noised_sigs = wf.read( '../data/noised.wav')noised_sigs = noised_sigs / 2 ** 15times = np.arange(len(noised_sigs)) / sample_ratefreqs = nf.fftfreq(times.size, 1 / sample_rate)noised_ffts = nf.fft(noised_sigs)noised_pows = np.abs(noised_ffts)fund_freq = np.abs(freqs[noised_pows.argmax()])print(fund_freq)noised_indices = np.where(np.abs(freqs) != fund_freq)filter_ffts = noised_ffts.copy()filter_ffts[noised_indices] = 0filter_pows = np.abs(filter_ffts)filter_sigs = nf.ifft(filter_ffts).realwf.write('../data/filter.wav', sample_rate, (filter_sigs * 2 ** 15).astype(np.int16))mp.figure('Filter', facecolor='lightgray')mp.subplot(221)mp.title('Time Domain', fontsize=16)mp.ylabel('Signal', fontsize=12)mp.tick_params(labelsize=10)mp.grid(linestyle=':')mp.plot(times[:178], noised_sigs[:178], c='orangered', label='Noised')mp.legend()mp.subplot(222)mp.title('Frequency Domain', fontsize=16)mp.ylabel('Power', fontsize=12)mp.tick_params(labelsize=10)mp.grid(linestyle=':')mp.semilogy(freqs[freqs >= 0], noised_pows[freqs >= 0], c='limegreen', label='Noised')mp.legend()mp.subplot(223)mp.xlabel('Time', fontsize=12)mp.ylabel('Signal', fontsize=12)mp.tick_params(labelsize=10)mp.grid(linestyle=':')mp.plot(times[:178], filter_sigs[:178], c='hotpink', label='Filter')mp.legend()mp.subplot(224)mp.xlabel('Frequency', fontsize=12)mp.ylabel('Power', fontsize=12)mp.tick_params(labelsize=10)mp.grid(linestyle=':')mp.plot(freqs[freqs >= 0], filter_pows[freqs >= 0], c='dodgerblue', label='Filter')mp.legend()mp.tight_layout()mp.show() 1000.0
随机数模块(random)
生成服从特定统计规律的随机数序列。
1. 二项分布
np.random.binomial(n, p, size)产生size个随机数，每个随机数来自n次尝试中的成功次数，其中每次尝试成功的概率为p。
# 猜硬币游戏：初始筹码1000，每轮猜9次，猜对5次及5次以上为赢，筹码加1，否则为输，筹码减1，求10000轮的过程中手中筹码的变化。# 代码：bi.pyimport numpy as npimport matplotlib.pyplot as mpoutcomes = np.random.binomial(9, 0.5, 10000)chips = [1000]for outcome in outcomes: if outcome >= 5: chips.append(chips[-1] + 1) else: chips.append(chips[-1] - 1)chips = np.array(chips)mp.figure('Binomial Distribution', facecolor='lightgray')mp.title('Binomial Distribution', fontsize=20)mp.xlabel('Round', fontsize=14)mp.ylabel('Chip', fontsize=14)mp.tick_params(labelsize=12)mp.grid(linestyle=':')o, h, l, c = 0, chips.argmax(), chips.argmin(), \ chips.size - 1if chips[o] < chips[c]: color = 'orangered'elif chips[c] < chips[o]: color = 'limegreen'else: color = 'dodgerblue'mp.plot(chips, c=color, label='Chip')mp.axhline(y=chips[o], linestyle='--', color='deepskyblue', linewidth=1)mp.axhline(y=chips[h], linestyle='--', color='crimson', linewidth=1)mp.axhline(y=chips[l], linestyle='--', color='seagreen', linewidth=1)mp.axhline(y=chips[c], linestyle='--', color='orange', linewidth=1)mp.legend()mp.show() 
2. 超几何分布
np.random.hypergeometric(ngood, nbad, nsample, size)产生size个随机数，每个随机数来自随机抽取nsample个样本中好样本的个数，总样本由ngood个好样本和nbad个坏样本组成
# 模球游戏：将25个好球和1个坏球放在一起，每次模3个球，全为好球加1分，只要摸到了坏球减6分，求100轮的过程中分值的变化。# 代码：hyper.pyimport numpy as npimport matplotlib.pyplot as mpoutcomes = np.random.hypergeometric(25, 1, 3, 100)scores = [0]for outcome in outcomes: if outcome == 3: scores.append(scores[-1] + 1) else: scores.append(scores[-1] - 6)scores = np.array(scores)mp.figure('Hypergeometric Distribution', facecolor='lightgray')mp.title('Hypergeometric Distribution', fontsize=20)mp.xlabel('Round', fontsize=14)mp.ylabel('Score', fontsize=14)mp.tick_params(labelsize=12)mp.grid(linestyle=':')o, h, l, c = 0, scores.argmax(), scores.argmin(), \ scores.size - 1if scores[o] < scores[c]: color = 'orangered'elif scores[c] < scores[o]: color = 'limegreen'else: color = 'dodgerblue'mp.plot(scores, c=color, label='Score')mp.axhline(y=scores[o], linestyle='--', color='deepskyblue', linewidth=1)mp.axhline(y=scores[h], linestyle='--', color='crimson', linewidth=1)mp.axhline(y=scores[l], linestyle='--', color='seagreen', linewidth=1)mp.axhline(y=scores[c], linestyle='--', color='orange', linewidth=1)mp.legend()mp.show() 
3. 标准正态分布
np.random.norm(size)`产生size个随机数，服从标准正态(平均值=0, 标准差=1)分布。
概率密度：$f(x) = \frac{
{ {e^{ - \frac{ { { {(x - \mu )}^2}}}{ {2{\sigma ^2}}}}}}}{ {\sigma \sqrt {2\pi } }}$
# norm.pyimport numpy as npimport matplotlib.pyplot as mpsamples = np.random.normal(size=10000)mp.figure('Normal Distribution', facecolor='lightgray')mp.title('Normal Distribution', fontsize=20)mp.xlabel('Sample', fontsize=14)mp.ylabel('Occurrence', fontsize=14)mp.tick_params(labelsize=12)mp.grid(axis='y', linestyle=':')bins = mp.hist(samples, 100, normed=True, edgecolor='steelblue', facecolor='deepskyblue', label='Normal')[1]probs = np.exp(-bins ** 2 / 2) / np.sqrt(2 * np.pi)mp.plot(bins, probs, 'o-', c='orangered', label='Probability')mp.legend()mp.show() 
杂项
1. 排序
联合间接排序
[张三 李四 王五 赵六] [70 60 80 70]<-[30 20 30 20]numpy.lexsort((参考序列, 待排序列))-> 有序索引numpy.sort_complex(复数数组)- 按照实部的升序排列，对于实部相同的元素，参考虚部的升序，直接返回排序后的结果数组。
numpy.searchsorted(有序序列, 待插序列)->位置数组- 表示将待插序列中的元素插入到有序序列中的哪些位置处，结果依然有序
numpy.insert(被插序列, 位置序列, 待插序列)->将待插序列中的元素- 按照位置序列中的位置，插入到被插序列中，返回插入后的结果
# sort.pyimport numpy as npages = np.array([30, 20, 30, 20])scores = np.array([70, 60, 80, 70])names = np.array(['zhangsan', 'lisi', 'wangwu', 'zhaoliu'])# 按照成绩的升序打印姓名，成绩相同的按照年龄的升序排列print(np.take(names, np.lexsort((ages, scores))))compleies = scores + ages * 1jprint(compleies)sorted_compleies = np.sort_complex(compleies)print(sorted_compleies)# 0 1 2 3 4 5 6a = np.array([1, 2, 4, 5, 6, 8, 9])b = np.array([7, 3])c = np.searchsorted(a, b)print(c)d = np.insert(a, c, b)print(d)
['lisi' 'zhaoliu' 'zhangsan' 'wangwu'][70.+30.j 60.+20.j 80.+30.j 70.+20.j][60.+20.j 70.+20.j 70.+30.j 80.+30.j][5 2][1 2 3 4 5 6 7 8 9]
2. 插值
import scipy.interpolate as sisi.interp1d(离散水平坐标, 离散垂直坐标, kind=插值算法(缺省为线性插值)) -> 插值器插值器(水平坐标)->垂直坐标
import numpy as npimport scipy.interpolate as siimport matplotlib.pyplot as mpmin_x, max_x = -2.5, 2.5con_x = np.linspace(min_x, max_x, 1001)con_y = np.sinc(con_x)dis_x = np.linspace(min_x, max_x, 11)dis_y = np.sinc(dis_x)# 线性插值linear = si.interp1d(dis_x, dis_y)lin_x = np.linspace(min_x, max_x, 51)lin_y = linear(lin_x)# 三次样条插值cubic = si.interp1d(dis_x, dis_y, kind='cubic')cub_x = np.linspace(min_x, max_x, 51)cub_y = cubic(cub_x)mp.figure('Interpolation', facecolor='lightgray')mp.subplot(221)mp.title('Continuous', fontsize=16)mp.ylabel('y', fontsize=12)mp.tick_params(labelsize=10)mp.grid(linestyle=':')mp.plot(con_x, con_y, c='hotpink', label='Continuous')mp.legend()mp.subplot(222)mp.title('Discrete', fontsize=16)mp.tick_params(labelsize=10)mp.grid(linestyle=':')mp.scatter(dis_x, dis_y, c='orangered', s=80, label='Discrete')mp.legend()mp.subplot(223)mp.title('Linear', fontsize=16)mp.xlabel('x', fontsize=12)mp.ylabel('y', fontsize=12)mp.tick_params(labelsize=10)mp.grid(linestyle=':')mp.plot(lin_x, lin_y, 'o-', c='limegreen', label='Linear')mp.scatter(dis_x, dis_y, c='orangered', s=80, zorder=3)mp.legend()mp.subplot(224)mp.title('Cubic', fontsize=16)mp.xlabel('x', fontsize=12)mp.tick_params(labelsize=10)mp.grid(linestyle=':')mp.plot(cub_x, cub_y, 'o-', c='dodgerblue', label='Cubic')mp.scatter(dis_x, dis_y, c='orangered', s=80, zorder=3)mp.legend()mp.tight_layout()mp.show() 
3. 积分
import scipy.integrate as sisi.quad(积分函数, 积分下限, 积分上限)->积分值, 最大误差
import numpy as npimport scipy.integrate as siimport matplotlib.pyplot as mpimport matplotlib.patches as mc# 定义函数def f(x): return 2 * x ** 2 + 3 * x + 4# 确定x范围a, b = -5, 5x1 = np.linspace(a, b, 1001)y1 = f(x1)# 用scipy来计算积分area = si.quad(f, a, b)[0]print(area)# 手动分割小梯形来计算面积，n越大，越接近真实值n = 10x2 = np.linspace(a, b, n + 1)y2 = f(x2)area = 0for i in range(n): area += (y2[i] + y2[i + 1]) * (x2[i + 1] - x2[i]) / 2print(area)# 绘制函数图mp.figure('Integral', dpi=120, facecolor='lightgray')mp.title('Integral', fontsize=20)mp.xlabel('x', fontsize=14)mp.ylabel('y', fontsize=14)mp.tick_params(labelsize=10)mp.grid(linestyle=':')mp.plot(x1, y1, c='orangered', linewidth=6, label=r'$y=2x^2+3x+4$', zorder=0)for i in range(n): mp.gca().add_patch(mc.Polygon([ [x2[i], 0], [x2[i], y2[i]], [x2[i + 1], y2[i + 1]], [x2[i + 1], 0]], fc='deepskyblue', ec='dodgerblue', alpha=0.5)) mp.legend()mp.show() 206.66666666666669210.0
4. 图像
scipy.ndimage中提供了一些简单的图像处理，如高斯模糊、任意角度旋转、边缘识别等功能。
import numpy as npimport scipy.misc as sm # 杂项工具import scipy.ndimage as snimport matplotlib.pyplot as mporiginal = sm.imread('../data/lily.jpg', True) # True就是黑白print(original.shape, original.dtype)median = sn.median_filter(original, (31, 31)) # 中值滤波，消除细节干扰、模糊化rotate = sn.rotate(original, 45)prewitt = sn.prewitt(original)mp.figure('Image', dpi=120, facecolor='lightgray')mp.subplot(221)mp.title('Original', fontsize=16)mp.axis('off')mp.imshow(original, cmap='gray')mp.subplot(222)mp.title('Median', fontsize=16)mp.axis('off')mp.imshow(median, cmap='gray')mp.subplot(223)mp.title('Rotate', fontsize=16)mp.axis('off')mp.imshow(rotate, cmap='gray')mp.subplot(224)mp.title('Prewitt', fontsize=16)mp.axis('off')mp.imshow(prewitt, cmap='gray')mp.tight_layout()mp.show() /Users/haoen110/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:6: DeprecationWarning: `imread` is deprecated!`imread` is deprecated in SciPy 1.0.0, and will be removed in 1.2.0.Use ``imageio.imread`` instead. (517, 690) float32
5. 金融
import numpy as np# 终值 = fv(利率, 期数, 每期支付, 现值)# 将1000元以1%的年利率存入银行5年，每年加存100元，# 到期后本息合计多少钱？fv = np.fv(0.01, 5, -100, -1000) # 凡是掏钱，都是负数print(round(fv, 2))# 现值 = pv(利率, 期数, 每期支付, 终值)# 将多少钱以1%的年利率存入银行5年，每年加存100元，# 到期后本息合计fv元？pv = np.pv(0.01, 5, -100, fv)print(pv)# 净现值 = npv(利率, 现金流)# 将1000元以1%的年利率存入银行5年，每年加存100元，# 相当于一次性存入多少钱？npv = np.npv(0.01, [-1000, -100, -100, -100, -100, -100])print(round(npv, 2))fv = np.fv(0.01, 5, 0, npv)print(round(fv, 2))# 内部收益率 = irr(现金流)# 将1000元存入银行5年，以后逐年提现100元、200元、# 300元、400元、500元，银行利率达到多少，可在最后# 一次提现后偿清全部本息，即净现值为0元？irr = np.irr([-1000, 100, 200, 300, 400, 500])print(round(irr, 2))npv = np.npv(irr, [-1000, 100, 200, 300, 400, 500])print(npv)# 每期支付 = pmt(利率, 期数, 现值)# 以1%的年利率从银行贷款1000元，分5年还清，# 平均每年还多少钱？pmt = np.pmt(0.01, 5, 1000)print(round(pmt, 2))# 期数 = nper(利率, 每期支付, 现值)# 以1%的年利率从银行贷款1000元，平均每年还pmt元，# 多少年还清？nper = np.nper(0.01, pmt, 1000)print(int(nper))# 利率 = rate(期数, 每期支付, 现值, 终值)# 从银行贷款1000元，平均每年还pmt元，nper年还清，# 年利率多少？rate = np.rate(nper, pmt, 1000, 0)print(round(rate, 2))
1561.11-1000.0-1485.341561.110.120.0-206.0450.01