python科学计算

一、NumPy-快速处理数据

1、NumPy提供了两种基本的对象

ndarray:存储单一数据类型的多维数组，称为数组

ufunc:对数组进行处理的特殊函数

2、查看numpy的版本

import numpy

numpy.version

3、创建数组

import numpy as np #引出np的对象
a=np.array([1,2,3,4])
b=np.array([5,6,7,8])
c=np.array([[1,2,3,4],[5,6,7,8],[7,8,9,10]])

a、shape获取数组的属性

>>> a.shape,b.shape,c.shape
((4,), (4,), (3, 4))
b、改变数组每个轴的长度
>>> c.shape=4,3
>>> c
array([[ 1,  2,  3],
[ 4,  5,  6],
[ 7,  8,  7],
[ 8,  9, 10]])
c、当某个元素的值为-1的时候，将自动计算此轴的长度
>>> c.shape=2,-1
>>> c
array([[ 1,  2,  3,  4,  5,  6],
[ 7,  8,  7,  8,  9, 10]])
d、reshape()创建指定形状的新数组
>>> d=a.reshape((2,2))#d=a.reshape(2,2)
>>> d
array([[1, 2],
[3, 4]])
>>> a
array([1, 2, 3, 4])
注：a、d同时共享数据存储空间
>>> a[1]=100
>>> a
array([  1, 100,   3,   4])
>>> d
array([[  1, 100],
[  3,   4]])

4、查看数据类型

c.dtype

dtype(‘int32’)

5、

二、Scipy——数值计算库

2.1常数和特殊函数

1、constants模块

from scipy import constants as C
print C.c#真空的光速
print C.h#普朗克常数
#字典physical_constants中，对应值分别为常数值、单位、误差
print C.physical_constants["electron mass"]#电子的质量
#单位信息
#1英里等于多少米，1英寸等于多少米，1克等于多少千克，1磅等于多少千克
print C.mile,C.inch,C.gram,C.pound

结果：

299792458.0

6.62606957e-34

(9.10938291e-31, 'kg', 4e-38)

1609.344 0.0254 0.001 0.45359237

2、special模块

伽马函数（gamma）：

import scipy.special as S
print S.gamma(4)
print S.gamma(0.5)
print S.gamma(1+1j)#支持复数
print S.gamma(1000)
print S.gammaln(1000)
结果：6.0
1.77245385091
(0.498015668118-0.154949828302j)
inf
5905.22042321

gamma函数阶乘和复数系上的扩展；gammaln是对其取对数

2.2拟合与优化——optimize

optimize提供了许多的数值优化算法

1、非线性方程组求解

fsolve()函数：基本调用形式：fsolve(func,x0)

func是计算方程组误差的函数；x时一个数组，其值为方程组的一组可能解。

方程组：f1(u1,u2,u3)=0,f2(u1,u2,u3)=0,f3(u1,u2,u3)=0,

那么func函数可以定义为：

def func(x):

u1,u2,u3=x

return [f1(u1,u2,u3),f2(u1,u2,u3),f3(u1,u2,u3)]

求解
5*x1+3=0,
4*x0*x0-2*sin(x1*x2)=0,
x1*x2-1.5=0
from math import sin
from scipy import optimize
def f(x):#计算方程组误差
x0,x1,x2=x.tolist()#调用数组的tolist方法，然后调用math的函数进行运算
#单个数值的计算中，标准浮点类型比Numpy浮点型要快许多，缩短时间
return [
5*x1+3,
4*x0*x0-2*sin(x1*x2),
x1*x2-1.5]
result=optimize.fsolve(f,[1,1,1])
#调用fsolve()时，传递计算误差的函数f()以及未知数的初始值
print result
print f(result)

雅可比矩阵：一阶偏导数以一定方式排列的矩阵



雅可比矩阵的函数j()和f()一样，类似

from math import sin,cos
from scipy import optimize
def f(x):
x0,x1,x2=x.tolist()
return [
5*x1+3,
4*x0*x0-2*sin(x1*x2),
x1*x2-1.5]
def  j(x):#计算方程组误差
x0,x1,x2=x.tolist()
return [
[0,5,0],
[8*x0,-2*x2*cos(x1*x2),-2*x1*cos(x1*x2)],
[0,x2,x1]
]
result = optimize.fsolve(f,[1,1,1],fprime=j)
#通过fprime=j参数将j()传递给fslove()
print result
print f(result)

2、最小二乘拟合：一组实验数据（x,y），通过参数p满足函数关系y=f(x),如f(x)=kx+b。

如果用p表示函数中需要确定的参数，则目标是找到一组p使得函数S的值最小:

import numpy as np
from scipy import optimize

X=np.array([8.19,2.72,6.39,8.71,4.7,2.66,3.78])
Y=np.array([7.01,2.78,6.47,6.71,4.1,4.23,4.05])
def residuals(p):
#计算以p为参数的直线和原始数据之间的误差
k,b = p
return Y-(k*X+ b)#注意缩进，否则出错
#leastsq使得residuals的输出数组的平方和最小，参数的初始值为【1，0】
r = optimize.leastsq(residuals,[1,0])
k,b=r[0]
print "k=",k,"b=",b

列二、`

import numpy as np

from scipy import optimize

import matplotlib.pyplot as plt

def logistic4(x, A, B, C, D):

return (A-D)/(1+(x/C)**B)+D

def residuals(p, y, x):

A, B, C, D = p

return y - logistic4(x, A, B, C, D)

def peval(x, p):

A, B, C, D = p

return logistic4(x, A, B, C, D)

A, B, C, D = .5, 2.5, 8, 7.3

x = np.linspace(0, 20, 20)

y_true = logistic4(x, A, B, C, D)

y_meas = y_true + 0.2 * np.random.randn(len(y_true))

p0 = [1/2]*4

plesq = optimize.leastsq(residuals, p0, args=(y_meas, x))

# leastsq函数的功能其实是根据误差（y_meas-y_true）

# 估计模型（也即函数）的参数

plt.figure(figsize=(6, 4.5))

plt.plot(x, peval(x, plesq[0]), x, y_meas, ‘o’, x, y_true)

plt.legend([‘Fit’, ‘Noisy’, ‘True’], loc=’upper left’)

plt.title(‘least square for the noisy data (measurements)’)

for i, (param, true, est) in enumerate(zip(‘ABCD’, [A, B, C, D], plesq[0])):

plt.text(11, 2-i*.5, ‘{} = {:.2f}, est({:.2f}) = {:.2f}’.format(param, true, param, est))

plt.savefig(‘./logisitic.png’)

plt.show()

`

列三：正弦波数据进行拟合

import numpy as np
#from math import sin
from scipy import optimize
import matplotlib.pyplot as pl
def func(x,p):
"""
数据拟合所用的函数:A*sin(2*pi*k*x+theta)
"""
A,k,theta=p
return A*np.sin(2*np.pi*k*x+theta)
def residuals(p,y,x):
"""
实验数据x,y和拟合函数之间的差，p为拟合需要找到的系数
"""
return y-func(x,p)
x=np.linspace(0,2*np.pi,100)
A,k,theta=10,0.34,np.pi/6#真实数据的函数参数
y0=func(x,[A,k,theta])#真实数据
#加入噪声之后的实验数据
np.random.seed(0)
y1=y0+2*np.random.randn(len(x))
p0=[7,0.40,0]#第一次猜测的函数拟合参数
#调用leastsq进行数据拟合
#residuals为计算误差的函数
#p0为拟合参数的初始值
#args为需要拟合的实验数据
plsq=optimize.leastsq(residuals,p0,args=(y1,x))
print u"真实参数：",[A,k,theta]
print u"拟合参数：",plsq[0]#实验数据拟合后的参数
pl.plot(x,y1,"o",label=u"带噪声的实验数据")
pl.plot(x,y0,label=u"真实数据")
pl.plot(x,func(x,plsq[0]),label=u"拟合数据")
pl.legend(loc="best")

带噪声的正弦波拟合