python做科学计算scipy

SciPy函数库在NumPy库的基础上增加了众多的数学、科学以及工程计算中常用的库函数。例如线性代数、常微分方程数值求解、信号处理、图像处理、稀疏矩阵等等。由于其涉及的领域众多、本书没有能力对其一一的进行介绍。作为入门介绍，让我们看看如何用SciPy进行插值处理、信号滤波以及用C语言加速计算。

3.1 最小二乘拟合

假设有一组实验数据(x[i], y[i])，我们知道它们之间的函数关系:y = f(x)，通过这些已知信息，需要确定函数中的一些参数项。例如，如果f是一个线型函数f(x) = k*x+b，那么参数k和b就是我们需要确定的值。如果将这些参数用 p 表示的话，那么我们就是要找到一组 p 值使得如下公式中的S函数最小：

这种算法被称之为最小二乘拟合(Least-square fitting)。

scipy中的子函数库optimize已经提供了实现最小二乘拟合算法的函数leastsq。下面是用leastsq进行数据拟合的一个例子：

# -*- coding: UTF-8 -*- 
'''
Created on 2017��8��24��
  
@author: zmz
'''
import numpy as np
from scipy.optimize import leastsq
from matplotlib import font_manager
import matplotlib.pyplot as plt;
  
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)


# 指定字体为宋体
zh_font = font_manager.FontProperties(fname=r'c:\windows\fonts\simsun.ttc', size=14)
x = np.linspace(0, -2*np.pi, 100)
A, k, theta = 10, 0.34, np.pi/6 # 真实数据的函数参数
y0 = func(x, [A, k, theta]) # 真实数据
y1 = y0 + 2 * np.random.randn(len(x)) # 加入噪声之后的实验数据    
  
p0 = [7, 0.2, 0] # 第一次猜测的函数拟合参数
  
# 调用leastsq进行数据拟合
# residuals为计算误差的函数
# p0为拟合参数的初始值
# args为需要拟合的实验数据
plsq = leastsq(residuals, p0, args=(y1, x))
  
print u"真实参数:", [A, k, theta] 
print u"拟合参数", plsq[0] # 实验数据拟合后的参数
  
plt.plot(x, y0, label=u"真实数据")
plt.plot(x, y1, label=u"带噪声的实验数据")
plt.plot(x, func(x, plsq[0]), label=u"拟合数据")
plt.legend(loc=(0.7, 0.8), prop=zh_font)
plt.show()

结果：

这个例子中我们要拟合的函数是一个正弦波函数，它有三个参数 A, k, theta ，分别对应振幅、频率、相角。假设我们的实验数据是一组包含噪声的数据 x, y1，其中y1是在真实数据y0的基础上加入噪声的到了。

通过leastsq函数对带噪声的实验数据x, y1进行数据拟合，可以找到x和真实数据y0之间的正弦关系的三个参数： A, k, theta。

真实参数: [10, 0.34, 0.5235987755982988]
拟合参数 [-10.15629453   0.34302702  -2.60935527]

我们看到拟合参数虽然和真实参数完全不同，但是由于正弦函数具有周期性，实际上拟合参数得到的函数和真实参数对应的函数是一致的。

3.2 函数最小值

optimize库提供了几个求函数最小值的算法：fmin, fmin_powell, fmin_cg, fmin_bfgs。

3.3 非线性方程组求解

optimize库中的fsolve函数可以用来对非线性方程组进行求解。

fsolve(func, x0)

func(x)是计算方程组误差的函数，它的参数x是一个矢量，表示方程组的各个未知数的一组可能解，func返回将x代入方程组之后得到的误差；x0为未知数矢量的初始值。如果要对如下方程组进行求解的话：

• 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 scipy.optimize import fsolve
from math import sin,cos


def f(x):
    x0 = float(x[0])
    x1 = float(x[1])
    x2 = float(x[2])
    return [
        5*x1+3,
        4*x0*x0 - 2*sin(x1*x2),
        x1*x2 - 1.5
    ]


result = fsolve(f, [1,1,1])


print result
print f(result)

结果：

[-0.70622057 -0.6        -2.5       ]
[0.0, -9.126033262418787e-14, 5.329070518200751e-15]

由于fsolve函数在调用函数f时，传递的参数为数组，因此如果直接使用数组中的元素计算的话，计算速度将会有所降低，因此这里先用float函数将数组中的元素转换为Python中的标准浮点数，然后调用标准math库中的函数进行运算。

在对方程组进行求解时，fsolve会自动计算方程组的雅可比矩阵，如果方程组中的未知数很多，而与每个方程有关的未知数较少时，即雅可比矩阵比较稀疏时，传递一个计算雅可比矩阵的函数将能大幅度提高运算速度。笔者在一个模拟计算的程序中需要大量求解近有50个未知数的非线性方程组的解。每个方程平均与6个未知数相关，通过传递雅可比矩阵的计算函数使计算速度提高了4倍。

3.7 滤波器设计

scipy.signal库提供了许多信号处理方面的函数。在这一节，让我们来看看如何利用signal库设计滤波器，查看滤波器的频率响应，以及如何使用滤波器对信号进行滤波。

假设如下导入signal库:

>>> import scipy.signal as signal

下面的程序设计一个带通IIR滤波器：

>>> b, a = signal.iirdesign([0.2, 0.5], [0.1, 0.6], 2, 40)

这个滤波器的通带为0.2*f0到0.5*f0，阻带为小于0.1*f0和大于0.6*f0，其中f0为1/2的信号取样频率，如果取样频率为8kHz的话，那么这个带通滤波器的通带为800Hz到2kHz。通带的最大增益衰减为2dB，阻带的最小增益衰减为40dB，即通带的增益浮动在2dB之内，阻带至少有40dB的衰减。

iirdesgin返回的两个数组b和a， 它们分别是IIR滤波器的分子和分母部分的系数。其中a[0]恒等于1。

3.8 用Weave嵌入C语言

Python作为动态语言其功能虽然强大，但是在数值计算方面有一个最大的缺点：速度不够快。在Python级别的循环和计算的速度只有C语言程序的百分之一。因此才有了NumPy, SciPy这样的函数库，将高度优化的C、Fortran的函数库进行包装，以供Python程序调用。如果这些高度优化的函数库无法实现我们的算法，必须从头开始写循环、计算的话，那么用Python来做显然是不合适的。因此SciPy提供了快速调用C++语言程序的方法-- Weave。下面是对NumPy的数组求和的例子：

解决方法;https://stackoverflow.com/questions/2817869/error-unable-to-find-vcvarsall-bat

和

http://blog.csdn.net/dlhlsc/article/details/62237595

重要：http://blog.csdn.net/secretx/article/details/17472107

问题描述：

Missing compiler_cxx fix for MSVCCompiler

distutils.errors.CompileError: error: Unable to find vcvarsall.bat

import scipy.weave as weave
import numpy as np
import time


def my_sum(a):
    n=int(len(a))
    code="""
    int i;


    double counter;
    counter =0;
    for(i=0;i<n;i++){
        counter=counter+a(i);
    }
    return_val=counter;
    """


    err=weave.inline(
        code,['a','n'],
        type_converters=weave.converters.blitz,
        compiler="gcc"
    )
    return err


a = np.arange(0, 10000000, 1.0)
# 先调用一次my_sum，weave会自动对C语言进行编译，此后直接运行编译之后的代码
my_sum(a)


start = time.clock()
for i in xrange(100):
    my_sum(a)  # 直接运行编译之后的代码
print "my_sum:", (time.clock() - start) / 100.0


start = time.clock()
for i in xrange(100):
    np.sum( a ) # numpy中的sum，其实现也是C语言级别
print "np.sum:", (time.clock() - start) / 100.0


start = time.clock()
print sum(a) # Python内部函数sum通过数组a的迭代接口访问其每个元素，因此速度很慢
print "sum:", time.clock() - start

结果：

在ubuntu上还是比较顺利的，

my_sum: 0.03872545
np.sum: 0.02853758
4.9999995e+13
sum: 2.994759

weave.inline函数的第一个参数为需要执行的C++语言代码，第二个参数是一个列表，它告诉weave要把Python中的两个变量a和n传递给C++程序，注意我们用字符串表示变量名。converters.blitz是一个类型转换器，将numpy的数组类型转换为C++的blitz类。C++程序中的变量a不是一个数组，而是blitz类的实例，因此它使用a(i)获得其各个元素的值，而不是用a[i]。最后我们通过compiler参数告诉weave要采用gcc为C++编译器。如果你安装的是python(x,y)的话，gcc(mingw32)也一起安装好了，否则你可能需要手工安装gcc编译器或者微软的Visual C++。

Note

 

在我的电脑上，虽然安装了Visual C++ 2008 Express，但仍然提示找不到合适的Visual C++编译器。似乎必须使用编译Python的编译器版本。因此还是用gcc来的方便。

本书的进阶部分还会对weave进行详细介绍。这段程序先给了我们一个定心丸：你再也不必担心Python的计算速度不够快了。