Python 科学计算方面
作者:Yupeng Jiang
- 伦敦大学学院 数学系 (英国顶尖大学,2018 QS世界大学排名中位列世界第7名,英国第3名)
- email:yupeng.jiang.13atcl.ac.uk
- 2016年6月5日
- [课件来自] https://zhuanlan.zhihu.com/p/21332075
遇到的模块
-
NumPy:多维数组的有效操作。 高效的数学函数。
-
Matplotlib:可视化:2D和(最近)3D图
-
SciPy:大型库实现各种数值算法,例如:
- 线性和非线性方程的解
- 优化
- 数值整合
-
Sympy:符号计算(解析的 Analytical)
-
Pandas:统计与数据分析(明天)
Numpy
Ndarray类型
NumPy提供了一种新的数据类型:ndarray(n维数组)。
- 与元组和列表不同,数组只能存储相同类型的对象
- 这使得数组上的操作比列表快得多; 此外,阵列占用的内存少于列表。
- 数组为列表索引机制提供强大的扩展。
均匀间隔的数组
import numpy as np # arange:range(start, stop, step)的所有三个参数 # 即起始值,结束值,步长都是可以用的 另外还有一个dtype参数,数据类型 a=np.arange(5) b=np.arange(10,100,20,dtype = float) #linspace(start,stop,num)返回数字间隔均匀的样本,按区间[start,stop]计算: c=np.linspace(0.,2.5,5)
多维数组矩阵
import numpy as np a = np.array([[1, 2, 3] ,[4, 5, 6]]) print(a) print(a.shape)#行,列数 print(a.ndim)#维度数 print(a.size)#元素个数
形状变化
import numpy as np a = np .arange(0, 20, 1) #1维 b = a.reshape((4, 5)) #4行5列 c = a.reshape((20, 1)) #2维 d = a.reshape((-1, 4)) #-1:自动确定 e = a.shape =(4, 5) #改变a的形状
Size(N,),(N,1)和(1,N)不同之处
import numpy as np a = np.array([1,2,3,4,5]) b = a.copy () c1 = np.dot(np.transpose(a), b)#矩阵乘法dot print(c1) c = np.array([[1,2,3],[4,5,6],[7,8,9]]) print(c) tran = c.transpose()#转置矩阵transpose print(tran) ax = np.reshape(a, (5,1)) bx = np.reshape(b, (1,5)) c = np.dot(ax, bx) print(c)
填充数组
import numpy as np a = np.zeros(3) b = np.zeros((2, 2), complex)#complex 复数 c = np.ones((2, 3)) d = np.random.rand(2, 4)#rand: 0和1之间均匀分布的随机数 e = np.random.randn(2, 4)#randn: 均值为0,标准差为1的标准(高斯)正态分布
副本和视图
- 采用标准的list切片建其副本
- 采用一个NumPy数组的切片可以在原始数组中创建一个视图。 两个数组都指向相同的内存。因此,当修改视图时,原始数组也被修改:
- 为避免修改原始数组,可以制作一个切片的副本
import numpy as np a = np.arange(5) print(a) b = a[2:];b[0] = 100 print(b) print(a) c = a[2:].copy(); c[0]=99 print(c) print(a)
矩阵乘法
- 运算符 * 表示元素乘法,而不是矩阵乘法:
- 使用dot()函数进行矩阵乘法:
- dot()方法也适用于矩阵向量(matrix-vector)乘法:
import numpy as np A = np.array([[1, 2],[3, 4]]) print(A * A) print(np.dot(A,A))
文件操作
- savetxt()将表保存到文本文件。
- 其他可用的格式(参见API文档)
- save()将表保存为Numpy“.npy”格式的二进制文件
- oadtxt()将以文本文件存储的表读入数组。
- 默认情况下,loadtxt()假定列是用空格分隔的。 您可以通过修改可选的参数进行更改。 以散列(#)开头的行将被忽略。
- 示例文本文件data.txt: (Example text file data.txt:) #Year Min temp. Hax temp.
- 1990 -1.5 25.3
- 1991 -3.2 21.2
import numpy as np a = np.linspace(0, 1, 12) a.shape=(3,4) np.savetxt('myfile.txt',a) np.save('myfile',a) table = np.loadtxt("data.txt")
Numpy包含更高效率的功能
-
Numpy包含许多常用的数学函数,例如:
- np.log
- np.maximum
- np.sin
- np.exp
- np.abs
-
在大多数情况下,Numpy函数比Math包中的类似函数更有效,特别是对于大规模数据。
Scipy
SciPy的结构
- scipy.integrate - >积分和普通微分方程
- scipy.linalg - >线性代数
- scipy.ndimage - >图像处理
- scipy.optimize - >优化和根查找(root finding)
- scipy.special - >特殊功能
- scipy.stats - >统计功能
- ...
要加载一个特定的模块,请这样使用, 例如 :
- from scipy import linalg
线性代数
线性方程的解
linalg的其他有用的方法:eig()(特征值和特征向量),det()(行列式)。
import numpy as np from scipy import linalg A = np.random.randn(5, 5) b = np.random.randn(5) x = linalg.solve(A, b) # A x = b #print(x) eigen = linalg.eig(A) # eigens #print(eigen) det = linalg.det(A) # determinant print(det)
数值整合
integration.quad是一维积分的自适应数值积分的函数。
import numpy as np from scipy import integrate def fun(x): return np.log(x) value, error = integrate.quad(fun,0,1) print(value) print(error)
用Scipy进行统计
scipy具有用于统计功能的子库,有一些有用的统计功能。 例如,给出标准正态分布的累积密度函数。
这个包,我们可以直接使用它,如下:
from scipy import stats y = stats.norm.cdf(1.2) #标准正态分布的累积密度函数 print(y)
优化:数据拟合
import numpy as np from scipy.optimize import curve_fit import matplotlib.pyplot as plt def func(x, a, b, c): return a * np.exp(-b * x) + c x = np.linspace(0, 4, 50) y = func(x, 2.5, 1.3, 0.5) ydata = y+0.2*np.random.normal(size=len(x)) popt, pcov = curve_fit(func, x, ydata) plt.plot(x, ydata, 'b*') plt.plot(x, func(x, popt[0],popt[1], popt[2]), 'r-') plt.title('$f(x)=ae^{-bx}+c$ curve fitting')
优化:根搜索
import numpy as np from scipy import optimize def fun(x): return np.exp(np.exp(x)) - x**2 # 通过初始化点0,找到兴趣0 (find zero of fun with initial point 0) # 通过Newton-Raphson方法 (by Newton-Raphson) value1 = optimize.newton(fun, 0) # 通过二分法找到介于(-5,5)之间的 (find zero between (-5,5) by bisection) value2 = optimize.bisect(fun, -5, 5)
Matplotlib
简单制图
import numpy as np import matplotlib.pyplot as plt x = np.linspace(0, 10, 201) #y = x ** 0.5 #plt.plot(x, y) # default plot plt.figure(figsize = (3, 3)) # new fig plt.plot(x, x**0.3, 'r--') # red dashed plt.plot(x, x-1, 'k-') # continue plot plt.plot(x, np.zeros_like(x), 'k-')
多个制图图例标签和标题
import numpy as np import matplotlib.pyplot as plt x = np.linspace(0, 10, 201) plt.figure(figsize = (4, 4)) for n in range(2, 5): y = x ** (1 / n) plt.plot(x, y, label='x^(1/' + str(n) + ')') plt.legend(loc = 'best') plt.xlabel('X axis') plt.ylabel('Y axis') plt.xlim(-2, 10) plt.title('Multi-plot e.g. ', fontsize = 18)
绘制子图
import numpy as np import matplotlib.pyplot as plt def pffcall(S, K): return np.maximum(S - K, 0.0) def pffput(S, K): return np.maximum(K - S, 0.0) S=np.linspace(50,151,100) fig = plt.figure(figsize=(12, 6)) sub1 = fig.add_subplot(121) sub1.set_title('Call', fontsize = 18) plt.plot(S, pffcall(S, 100), 'r-', lw = 4)#lw是宽度 plt.plot(S, np.zeros_like(S), 'black',lw = 1) sub1.grid(True)#显示网格线 sub1.set_xlim([60, 120]) sub1.set_ylim([-10, 40]) sub2 = fig.add_subplot(122) sub2.set_title('Put', fontsize = 18) plt.plot(S, pffput(S, 100), 'r-', lw = 4) plt.plot(S, np.zeros_like(S), 'black',lw = 1) sub2.grid(True) sub2.set_xlim([60, 120]) sub2.set_ylim([-10, 40])
在绘制的图上添加文本和注释
import numpy as np from scipy.stats import norm import matplotlib.pyplot as plt def call(S, K=100, T=0.5, vol=0.6, r=0.05): d1 = (np.log(S/K) + (r + 0.5 * vol**2) \ *T) / np.sqrt(T) / vol d2 = (np.log(S/K) + (r - 0.5 * vol**2) \ *T) / np.sqrt(T) / vol return S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2) def delta(S, K=100, T=0.5, vol=0.6, r=0.05): d1 = (np.log(S/K) + (r + 0.5 * vol**2)*T) / np.sqrt(T) / vol return norm.cdf(d1) S = np.linspace(40, 161, 100) fig = plt.figure(figsize=(7, 6)) ax = fig.add_subplot(111) #add_subplot(111)把画布分为1行1列,area为从左往右从上往下第1块区域 plt.plot(S,(call(S)-call(100)),'r',lw=1) plt.plot(100, 0, 'ro', lw=1) plt.plot(S,np.zeros_like(S), 'black', lw = 1) plt.plot(S,call(S)-delta(100)*S- \ (call(100)-delta(100)*100), 'y', lw = 1) ax.annotate('$\Delta$ hedge', xy=(100, 0), \ xytext=(110, -10),arrowprops= \ dict(headwidth =3,width = 0.5, \ facecolor='black', shrink=0.05))#xy是点的位置 #xytext是文字位置,arrowprops 箭头 ax.annotate('Original call', xy= \ (120,call(120)-call(100)),xytext\ =(130,call(120)-call(100)),\ arrowprops=dict(headwidth =10,\ width = 3, facecolor='cyan', \ shrink=0.05)) plt.grid(True) plt.xlim(40, 160) plt.xlabel('Stock price', fontsize = 18) plt.ylabel('Profits', fontsize = 18)
3D制图
import numpy as np import matplotlib.pyplot as plt from matplotlib import cm from mpl_toolkits.mplot3d import Axes3D x, y = np.mgrid[-5:5:100j, -5:5:100j] z = x**2 + y**2 fig = plt.figure(figsize=(8, 6)) ax = plt.axes(projection='3d') surf = ax.plot_surface(x, y, z, rstride=1,\ cmap=cm.coolwarm, cstride=1, \ linewidth=0) fig.colorbar(surf, shrink=0.5, aspect=5) plt.title('3D plot of $z = x^2 + y^2$')
3D制图练习
import numpy as np import matplotlib.pyplot as plt from matplotlib import cm from mpl_toolkits.mplot3d import Axes3D x, y = np.mgrid[-10:10:100j, -10:10:100j] z = x**2*0.05 + np.sin(x) + y**2*0.05 + np.sin(y) fig = plt.figure(figsize=(8, 6)) ax = plt.axes(projection='3d') surf = ax.plot_surface(x, y, z, rstride=1,\ cmap=cm.rainbow, cstride=1, \ linewidth=0) fig.colorbar(surf, shrink=0.5, aspect=5) plt.title('3D plot of $z = 0.05*x^2 + sin(x) + 0.05*y^2+sin(y)$')
♥形线
import numpy as np import matplotlib.pyplot as plt def y(S): return 13*np.cos(S)-5*np.cos(2*S)-2*np.cos(3*S)-np.cos(4*S) def x(S): return 16*np.sin(S)**3 S = np.linspace(0, 2*np.pi, 100) fig = plt.figure(figsize=(4, 4)) plt.plot(x(S),y(S),'blue',lw=1) plt.title("Heart") plt.grid(True)
符号计算(Symbolic computat)
以上只是数值计算,接下来涉及符号计算,python通过模块sysmpy来进行符号计算,类似于方程求解,积分等的显式求解。
声明一个符号变量
import sympy as sy #声明x,y为变量 x = sy.Symbol('x') y = sy.Symbol('y') a, b = sy.symbols('a b') #创建一个新符号(不是函数 f = x**2 + y**2 -2*x*y + 5 print(f) #自动简化 g = x**2 + 2 - 2*x + x**2 -1 print(g)
符号的使用1:求解方程
import sympy as sy x = sy.Symbol ('x') y = sy.Symbol('y') # 给定[-1,1] (give [-1, 1]) print(sy.solve (x**2 - 1)) # 无解 (no guarantee for solution) print(sy.solve(x**3 + 0.5*x**2 - 1)) # 用x的表达式表示y (exepress x in terms of y) print (sy.solve(x**3 + y**2)) # 错误:找不到算法 (error: no algorithm can be found) print(sy.solve(x**x + 2*x - 1))
符号的使用2:集成
import sympy as sy x = sy.Symbol('x') y = sy.Symbol( 'y') a,b = sy.symbols ( 'a b') # 单变量 single variable f = sy.sin(x) + sy.exp(x) print(sy.integrate(f, (x, a, b))) print(sy.integrate(f, (x, 1, 2))) print(sy.integrate(f, (x, 1.0,2.0))) # 多变量 multi variables g = sy.exp(x) + x * sy.sin(y) print(sy.integrate(g, (y,a,b)))
符号的使用3:分化
import sympy as sy x = sy.Symbol( 'x') y = sy.Symbol( 'y') # 单变量 (single variable) f = sy.cos(x) + x**x print(sy . diff (f , x)) # 多变量 (multi variables) g = sy.cos(y) * x + sy.log(y) print(sy.diff (g, y))
代码整合
链接:https://pan.baidu.com/s/1X6eG5C03KmBJ8yN-WoRRSw
提取码:a6dy
复制这段内容后打开百度网盘手机App,操作更方便哦