3.Python scipy笔记（进行科学计算）

2018-12-04  本文已影响0人  Louis_Duan

1 简介

1.1 定义

scipy是一个以numpy为基础的用来进行科学计算的一个第三方库。

1.2 构成

子模块           描述
cluster         聚类算法
constants       物理数学常数
fftpack         快速傅里叶变换
integrate       积分和常微分方程求解
interpolate     插值
io              输入输出
linalg          线性代数
odr             正交距离回归
optimize        优化和求根
signal          信号处理
sparse          稀疏矩阵
spatial         空间数据结构和算法
special         特殊方程
stats           统计分布和函数
weave           C/C++ 积分

2 插值

2.1 一维插值

from scipy.interpolate import interp1d
interp1d(x,y,kind='linear',...)

参数 x,y是一系列已知的点，kind参数是插值的类型

候选值 作用
‘zero’、‘nearest’ 0阶插值、最近邻插值
‘slinear’、‘linear’ 线性插值，用一条直线连接所有取样点，相当于一阶B样条曲线
‘quadratic’、‘cubic’ 二次插值，更高阶的曲线可以直接使用整数值指定

2.2 代码示例

import numpy as np
from scipy.interpolate import interp1d

#创建待插值的数据
x = np.linspace(0, 10*np.pi, 20)
y = np.cos(x)

# 分别用linear和quadratic插值
fl = interp1d(x, y, kind='linear')
fq = interp1d(x, y, kind='quadratic')

#设置x的最大值和最小值以防止插值数据越界
xint = np.linspace(x.min(), x.max(), 1000)
yintl = fl(xint)
yintq = fq(xint)

作者：火锅侠
链接：https://www.jianshu.com/p/b306095309db
來源：简书
简书著作权归作者所有，任何形式的转载都请联系作者获得授权并注明出处。

2.3 径向基函数

径向基函数，简单来说就是 x 处的函数值只依赖于 x 和 某点 c 的距离。

2.3.1 常用的径向基（RBF）函数有：

高斯函数

x = np.linspace(-3,3,100)
plt.plot(x,np.exp(-1*x**2))
plt.title(''Gaussian'')
Figure_1.png

Mutiquadric函数（多元二次函数）

plt.plot(x,np.sqrt(1+x**2))
plt.title(''Multiquadric'')
Figure_1.png

Inverse quadratic函数（逆二次函数）

plt.plot(x,1./np.sqrt(1+x**2))
plt.title(''Inverse Multiquadric'')
Figure_1.png

2.3.2 径向基函数插值

from scipy.interpolate.rbf import Rbf
# 使用multiquadric的核
cp_rbf = Rbf(x,y,function="multiquadric")
# 使用gaussian的核
cp_rbf=Rbf(x,y,function="gaussian")
# 使用inverse_multiquadric的核
cp_rbf=Rbf(x,y,function="inverse_multiquadric")

2.3.3 高维RBF插值

import numpy as np
import pylab as pl
from  mpl_toolkits.mplot3d import Axes3D
from scipy.interpolate.rbf import Rbf

# 3 维数据点
x,y=np.mgrid[-np.pi/2:np.pi/2:5j, -np.pi/2:np.pi/2:5j]
z = np.cos(np.sqrt(x**2+y**2))
fig = pl.figure(figsize=(12,6))
ax = fig.gca(projection="3d")
ax.scatter(x,y,z)

# 3 维RBF插值
zz = Rbf(x,y,z)
xx, yy = np.mgrid[-np.pi/2:np.pi/2:50j, -np.pi/2:np.pi/2:50j]
fig = pl.figure(figsize=(12,6))
ax = fig.gca(projection="3d")
ax.plot_surface(xx,yy,zz(xx,yy),rstride=1,cstride=1,cmap=pl.cm.jet)

pl.show()

3 概率统计方法

4 曲线拟合

什么是曲线拟合：找一个方程来拟合数据。
基础包

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from numpy import polyfit,poly1d
# polyfit返回的是系数，poly1d返回的函数式

4.1 多项式拟合

从一阶到九阶拟合多项式拟合正弦函数

import numpy as np
import matplotlib.pyplot as plt
from numpy import  polyfit,poly1d
x = np.linspace(-np.pi,np.pi,100)
y = np.sin(x)
y1 = poly1d(polyfit(x,y,1))
y3 = poly1d(polyfit(x,y,3))
y5 = poly1d(polyfit(x,y,5))
y7 = poly1d(polyfit(x,y,7))
y9 = poly1d(polyfit(x,y,9))
x = np.linspace(-3 * np.pi,3 * np.pi,100)

p = plt.plot(x, np.sin(x), 'k')
p = plt.plot(x, y1(x),label='y1')
p = plt.plot(x, y3(x),label='y3')
p = plt.plot(x, y5(x),label='y5')
p = plt.plot(x, y7(x),label='y7')
p = plt.plot(x, y9(x),label='y9')

a = plt.axis([-3 * np.pi, 3 * np.pi, -1.25, 1.25])

plt.legend()
plt.show()
Figure_1.png

4.2 最小二乘拟合

什么是最小二乘：已知多个近似分布于同一直线上的点，想拟合出一个直线方程：设该直线方程为y=kx+b，调整参数k和b，使得所有点到该直线的距离平方之和最小，设此时满足要求的k=k0，b=b0，则直线方程为y=k0x+b0。

相关模块

from scipy.linalg import lstsq#求最小二乘解
from scipy.stats import linregress#线性回归
from scipy.optimize import leastsq
#（func，x0，args=（））
func是自己定义的一个计算误差的函数
x0是计算的初始参数值
args是指定func的其他参数

用leastsq函数拟合数据，待拟合函数为kx+b

###最小二乘法试验###
import numpy as np
from scipy.optimize import leastsq

###采样点(Xi,Yi)###
Xi=np.array([8.19,2.72,6.39,8.71,4.7,2.66,3.78])
Yi=np.array([7.01,2.78,6.47,6.71,4.1,4.23,4.05])

###需要拟合的函数func及误差error###
def func(p,x):
    k,b=p
    return k*x+b

def error(p,x,y,s):
    print (s)
    return func(p,x)-y #x、y都是列表，故返回值也是个列表

#TEST
p0=[100,2]
#print( error(p0,Xi,Yi) )

###主函数从此开始###
s="Test the number of iteration" #试验最小二乘法函数leastsq得调用几次error函数才能找到使得均方误差之和最小的k、b
Para=leastsq(error,p0,args=(Xi,Yi,s)) #把error函数中除了p以外的参数打包到args中
k,b=Para[0]
print("k=",k,'\n',"b=",b)

###绘图，看拟合效果###
import matplotlib.pyplot as plt

plt.figure(figsize=(8,6))
plt.scatter(Xi,Yi,color="red",label="Sample Point",linewidth=3) #画样本点
x=np.linspace(0,10,1000)
y=k*x+b
plt.plot(x,y,color="orange",label="Fitting Line",linewidth=2) #画拟合直线
plt.legend()
plt.show()
Figure_1.png

5 最小化函数

？？

6 积分

？？

7 解微分方程

？？

8 稀疏矩阵

Scipy 提供了稀疏矩阵的支持（scipy.sparse）。稀疏矩阵主要使用 位置 + 值 的方法来存储矩阵的非零元素
常用的存储格式有：

from scipy.sparse import *
import numpy as np
A = coo_matrix([[1,2,0],[0,0,3],[4,0,5]])
print (A)
  (0, 0)    1
  (0, 1)    2
  (1, 2)    3
  (2, 0)    4
  (2, 2)    5
前面是非零值的坐标，后面是值
A.toarray()
array([[1, 2, 0],
       [0, 0, 3],
       [4, 0, 5]])

C 
  (0, 0)    3
  (0, 2)    1
  (1, 1)    2
  (3, 3)    1
from scipy import sparse

row, col, val = sparse.find(C)
print (row, col, val)
[0 0 1 3] [0 2 1 3] [3 1 2 1]

9 线性代数

numpy 和 scipy 中，负责进行线性代数部分计算的模块叫做 linalg。

线性代数的基本操作对象是矩阵

numpy.matrix创建矩阵

A = np.mat("[1, 2; 3, 4]")
print (repr(A))
matrix([[1, 2],
        [3, 4]])

转置和逆

A.I
A.T

矩阵乘法

A*B

2维numpy.ndarray创建矩阵

A = np.array([[1,2], [3,4]])
>> array([[1, 2],
       [3, 4]])

逆

linalg.inv(A)

转置

A.T

矩阵乘法

A.dot(b)

普通乘法

A*B

求解线性方程组

x+3y+5z=10
2x+5y+z=8
2x+3y+8z=3

>>A = np.array([[1, 3, 5],
              [2, 5, 1],
              [2, 3, 8]])
>>b = np.array([10, 8, 3])
>>linalg.solve(A,b)
[-9.28  5.16  0.76]

计算行列式

>>A = np.array([[1, 3, 5],
              [2, 5, 1],
              [2, 3, 8]])
>>linalg.det(A)
>>-25.0

计算矩阵或向量的模

x_norm=np.linalg.norm(x, ord=None, axis=None, keepdims=False)

向量范数：

参数 说明 计算方法
ord=2 二范数 每个元素的平方的和再开方
ord=1 一范数 每个元素的绝对值的和
ord=np.inf 无穷范数 max（|xi|） 
>>A = np.array([[1, 2],
              [3, 4]])
>>linalg.norm(A)
import numpy as np
x = np.array([
    [0, 3, 4],
    [1, 6, 4]])
print linalg.norm(A)

print linalg.norm(A,'fro') # frobenius norm 默认值

print linalg.norm(A,1) # L1 norm 最大列和

print linalg.norm(A,-1) # L -1 norm 最小列和

print linalg.norm(A,np.inf) # L inf norm 最大行和

伪逆

linalg.pinv:使用最小二乘的解法求伪逆
linalg.pinv2:使用求奇异值的算法求解

特征值和特征向量

linalg.eig(A):返回矩阵的特征值和特征向量
linalg.eigvals(A):返回矩阵的特征值
linalg.eig(A,B):?  +
上一篇下一篇

猜你喜欢

热点阅读