3.Python scipy笔记(进行科学计算)

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)。稀疏矩阵主要使用 位置 + 值 的方法来存储矩阵的非零元素
常用的存储格式有:

  • COO 格式在构建矩阵时比较高效
  • CSC 和 CSR 格式在乘法计算时比较高效
    例子:
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]])
  • sparse.find()
    返回一个三元组,表示稀疏矩阵中的非零元素的(行,列,值)
    例子:
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]
  • sparse.issparse()
    查看一个矩阵是否为稀疏矩阵,返回值为True或者False

9 线性代数

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

  • numpy.linalg(推荐使用)和scipy.linalg 区别:
    1.scipy.linalg包含scipy.linalg中所有函数,同时还包含了很多其没有的函数。
    2.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|)
  • 矩阵的范数
    ord=1:列和范数,即所有矩阵列向量绝对值之和的最大值
    ord=2:谱范数,即ATA矩阵的最大特征值的开平方。
    ord=inf:行和范数,即所有矩阵行向量绝对值之和的最大值.
  • axis:处理类型
    axis=1表示按行向量处理,求多个行向量的范数
    axis=0表示按列向量处理,求多个列向量的范数
    axis=None表示矩阵范数。
>>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):?  +
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 203,324评论 5 476
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 85,303评论 2 381
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 150,192评论 0 337
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,555评论 1 273
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,569评论 5 365
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,566评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,927评论 3 395
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,583评论 0 257
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,827评论 1 297
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,590评论 2 320
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,669评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,365评论 4 318
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,941评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,928评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,159评论 1 259
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,880评论 2 349
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,399评论 2 342

推荐阅读更多精彩内容

  • 原文章为scikit-learn中"用户指南"-->"监督学习的第七节:Gaussian Processes"##...
    HabileBadger阅读 18,295评论 0 9
  • 梦里有海的,梦里真的有海,一望无际的大海,看不到尽头,我猜应该是差不多傍晚了。 我们并排光着脚丫站在还有一丝暖热的...
    susumini阅读 791评论 0 0
  • Spoon.bat:图形界面方式启动作业和转换设计器 Pan.bat:命令行方式执行转换 Kitchen.bat:...
    laddie_592a阅读 1,333评论 0 2
  • 生活中我们总会遇到一些不会聊天的人。你跟他说“这个才好吃吗”,他说“还行”这样本来想说话的你瞬间就没有什么话题想说...
    微风知意阅读 237评论 0 1