数学建模三剑客 MSN
(点击上方公号，快速关注我们)
作者：许向武
前言
不管是不是巴萨的球迷，只要你喜欢足球，就一定听说过梅西(Messi)、苏亚雷斯(Suarez)和内马尔(Neymar)这个MSN组合。在众多的数学建模辅助工具中，也有一个犀利无比的MSN组合，他们就是python麾下大名鼎鼎的 Matplotlib + Scipy + Numpy三剑客。
本文是我整理的MSN学习笔记，有些理解可能比较肤浅，甚至是错误的。如果因此误导了某位看官，在工作中造成重大失误或损失，我顶多只能赔偿一顿饭——还得是我们楼下的十元盒饭。特此声明。
文中代码均从我的这台时不时出点问题、闹个情绪的Yoga 3 pro上复制而来，这意味着所有的代码均可在下面的运行环境中顺利运行：
pyhton 2.7.8
numpy 1.11.1
scipy 0.16.1
matplotlib 1.5.1
三剑客之Numpy
numpy是一个开源的python科学计算库，包含了很多实用的数学函数，涵盖线性代数、傅里叶变换和随机数生成等功能。最初的numpy其实是scipy的一部分，后来才从scipy中分离出来。
numpy不是python的标准库，需要单独安装。假定你的运行环境已经安装了python包管理工具pip，numpy的安装就非常简单：
pip install numpy
数组对象
ndarray是多维数组对象，也是numpy最核心的对象。在numpy中，数组的维度(dimensions)叫做轴(axes)，轴的个数叫做秩(rank)。通常，一个numpy数组的所有元素都是同一种类型的数据，而这些数据的存储和数组的形式无关。
下面的例子，创建了一个三维的数组（在导入numpy时，一般都简写成np）。
import numpy as np
a = np.array([[1,2,3],[4,5,6],[7,8,9]])
数据类型
numpy支持的数据类型主要有布尔型(bool)、整型(integrate)、浮点型(float)和复数型(complex)，每一种数据类型根据占用内存的字节数又分为多个不同的子类型。常见的数据类型见下表。
创建数组
通常，我们用np.array()创建数组。如果仅仅是创建一维数组，也可以使用np.arange()或者np.linspace()的方法。np.zeros()、np.ones()、np.eye()则可以构造特殊的数据。np.random.randint()和np.random.random()则可以构造随机数数组。
构造复杂数组
很多时候，我们需要从简单的数据结构，构造出复杂的数组。例如，用一维的数据生成二维格点。
重复数组: tile
一维数组网格化: meshgrid
指定范围和分割方式的网格化: mgrid
上面的例子中用到了虚数。构造虚数的方法如下：
>>> complex(2,5)
(2+5j)
数组的属性
numpy的数组对象除了一些常规的属性外，也有几个类似转置、扁平迭代器等看起来更像是方法的属性。扁平迭代器也许是遍历多维数组的一个简明方法，下面的代码给出了一个例子。
改变数组维度
numpy数组的存储顺序和数组的维度是不相干的，因此改变数组的维度是非常便捷的操作，除resize()外，这一类操作不会改变所操作的数组本身的存储顺序。
索引和切片
对于一维数组的索引和切片，numpy和python的list一样，甚至更灵活。
假设有一栋2层楼，每层楼内的房间都是3排4列，那我们可以用一个三维数组来保存每个房间的居住人数（当然，也可以是房间面积等其他数值信息）。
数组合并
数组合并除了下面介绍的水平合并、垂直合并、深度合并外，还有行合并、列合并，以及concatenate()等方式。假如你比我还懒，那就只了解前三种方法吧，足够用了。
数组拆分
拆分是合并的逆过程，概念是一样的，但稍微有一点不同：
数组运算
数组和常数的四则运算，是数组的每一个元素分别和常数运算；数组和数组的四则运算则是两个数组对应元素的运算（两个数组有相同的shape，否则抛出异常）。
特别提示：如果想对数组内符合特定条件的元素做特殊处理，下面的代码也许有用。
数组方法和常用函数
数组对象本身提供了计算算数平均值、求最大最小值等内置方法，numpy也提供了很多实用的函数。为了缩减篇幅，下面的代码仅以一维数组为例，展示了这些方法和函数用法。事实上，大多数情况下这些方法和函数对于多维数组同样有效，只有少数例外，比如compress函数。
矩阵对象
matrix是矩阵对象，继承自ndarray类型，因此含有ndarray的所有数据属性和方法。不过，当你把矩阵对象当数组操作时，需要注意以下几点：
matrix对象总是二维的，即使是展平（ravel函数）操作或是成员选择，返回值也是二维的
matrix对象和ndarray对象混合的运算总是返回matrix对象
创建矩阵
matrix对象可以使用一个Matlab风格的字符串来创建（以空格分隔列，以分号分隔行的字符串），也可以用数组来创建。
矩阵的特有属性
矩阵有几个特有的属性使得计算更加容易，这些属性有：
矩阵乘法
对ndarray对象而言，星号是按元素相乘，dot()函数则当作矩阵相乘。对于matrix对象来说，星号和dot()函数都是矩阵相乘。特别的，对于一维数组，dot()函数实现的是向量点乘（结果是标量），但星号实现的却不是差乘。
线性代数模块
numpy.linalg 是numpy的线性代数模块，可以用来解决逆矩阵、特征值、线性方程组以及行列式等问题。
计算逆矩阵
尽管matrix对象本身有逆矩阵的属性，但用numpy.linalg模块求解矩阵的逆，也是非常简单的。
计算行列式
如何计算行列式，我早已经不记得了，但手工计算行列式的痛苦，我依然记忆犹新。现在好了，你在手机上都可以用numpy轻松搞定（前提是你的手机上安装了python + numpy）。
m = np.mat('0 1 2; 1 0 3; 4 -3 8')
np.linalg.det(m) # 什么？这就成了？
2.0
计算特征值和特征向量
截至目前，我的工作和特征值、特征向量还有没任何关联。记录这一节，纯粹是为了我女儿，她正在读数学专业。
求解线性方程组
有线性方程组如下：
x – 2y + z = 0
2y -8z = 8
-4x + 5y + 9z = -9
求解过程如下：
三剑客之Matplotlib
matplotlib 是python最著名的绘图库，它提供了一整套和Matlab相似的命令API，十分适合交互式地进行制图。而且也可以方便地将它作为绘图控件，嵌入GUI应用程序中。matplotlib 可以绘制多种形式的图形包括普通的线图，直方图，饼图，散点图以及误差线图等；可以比较方便的定制图形的各种属性比如图线的类型，颜色，粗细，字体的大小等；它能够很好地支持一部分 TeX 排版命令，可以比较美观地显示图形中的数学公式。
pyplot介绍
Matplotlib 包含了几十个不同的模块， 如 matlab、mathtext、finance、dates 等，而 pyplot 则是我们最常用的绘图模块，这也是本文介绍的重点。
中文显示问题的解决方案
有很多方法可以解决此问题，但下面的方法恐怕是最简单的解决方案了（我只在windows平台上测试过，其他平台请看官自测）。
绘制最简单的图形
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> x = np.arange(0, 2*np.pi, 0.01)
>>> y = np.sin(x)
>>> plt.plot(x, y)
>>> plt.show()
设置标题、坐标轴名称、坐标轴范围
如果你在python的shell中运行下面的代码，而shell的默认编码又不是utf-8的话，中文可能仍然会显示为乱码。你可以尝试着把 u’正弦曲线’ 写成 ‘正弦曲线’.decode(‘gbk’)或者‘正弦曲线’.decode(‘utf-8’)
设置点和线的样式、宽度、颜色
plt.plot函数的调用形式如下：
plot(x, y, color='green', linestyle='dashed', linewidth=1, marker='o', markerfacecolor='blue',markersize=6)
plot(x, y, c='g', ls='--', lw=1, marker='o', mfc='blue', ms=6)
color指定线的颜色，可简写为“c”。颜色的选项为：
蓝色： ‘b’ (blue)
绿色： ‘g’ (green)
红色： ‘r’ (red)
墨绿： ‘c’ (cyan)
洋红： ‘m’ (magenta)
黄色： ‘y’ (yellow)
黑色： ‘k’ (black)
白色： ‘w’ (white)
灰度表示： e.g. 0.75 ([0,1]内任意浮点数)
RGB表示法： e.g. ‘#2F4F4F’ 或 (0.18, 0.31, 0.31)
linestyle指定线型，可简写为“ls”。线型的选项为：
实线： ‘-’ (solid line)
虚线： ‘–’ (dashed line)
虚点线： ‘-.’ (dash-dot line)
点线： ‘:’ (dotted line)
无： ”或’ ‘或’None’
linewidth指定线宽，可简写为“lw”。
marker描述数据点的形状
点线： ‘.’
点线： ‘o’
加号： ‘+
叉号： ‘x’
上三角： ‘^’
上三角： ‘v’
markerfacecolor指定数据点标记的表面颜色，可 简写为“ mfc”。
markersize指定数据点标记的大小，可 简写为“ ms”。
文本标注和图例
我们分别使用不同的线型、颜色来绘制以10、e、2为基的一组幂函数曲线，演示文本标注和图例的使用。
在绘制图例时，loc用于指定图例的位置，可用的选项有：
best
upper right
upper left
lower left
lower right
绘制多轴图
在介绍如何将多幅子图绘制在同一画板的同时，顺便演示如何绘制直线和矩形。我们可以使用subplot函数快速绘制有多个轴的图表。subplot函数的调用形式如下：
subplot(numRows, numCols, plotNum)
subplot将整个绘图区域等分为numRows行 * numCols列个子区域，然后按照从左到右，从上到下的顺序对每个子区域进行编号，左上的子区域的编号为1。如果numRows，numCols和plotNum这三个数都小于10的话，可以把它们缩写为一个整数，例如subplot(323)和subplot(3,2,3)是相同的。subplot在plotNum指定的区域中创建一个轴对象。如果新创建的轴和之前创建的轴重叠的话，之前的轴将被删除。
三剑客之Scipy
前面已经说过，最初的numpy其实是scipy的一部分，后来才从scipy中分离出来。scipy函数库在numpy库的基础上增加了众多的数学、科学以及工程计算中常用的库函数。例如线性代数、常微分方程数值求解、信号处理、图像处理、稀疏矩阵等等。由于其涉及的领域众多，我之于scipy，就像盲人摸大象，只能是摸到哪儿算哪儿。
插值
一维插值和二维插值，是我最常用的scipy的功能之一，也是最容易上手的。
一维插值和样条插值
下面的例子清楚地展示了线性插值和样条插值之后的数据形态。
将原始数据以及线性插值和样条插值之后的数据绘制在一起，效果会比较明显：
代码如下：
特别说明：样条插值附带了很多默认参数，下面是简单的说明。详情请自行搜索。
scipy.interpolate.splrep(x, y, w=None, xb=None, xe=None, k=3, task=0, s=None, t=None, full_output=0,per=0, quiet=1)
# 参数s用来确定平滑点数，通常是m-SQRT(2m),m是曲线点数。如果在插值中不需要平滑应该设定s=0。splrep()输出的是一个3元素的元胞数组（t,c,k）,其中t是曲线点，c是计算出来的系数，k是样条阶数，通常是3阶，但可以通过k改变。
scipy.interpolate.splev(x, tck, der=0)
# 其中的der是进行样条计算是需要实际计算到的阶数，必须满足条件der<=k
拟合
在工作中，我们常常需要在图中描绘某些实际数据观察的同时，使用一个曲线来拟合这些实际数据。所谓拟合，就是找出符合数据变化趋势的曲线方程，或者直接绘制出拟合曲线。
使用numpy.polyfit拟合
下面这段代码，基于Numpy模块，可以直接绘制出拟合曲线，但我无法得到曲线方程（尽管输出了一堆曲线参数）。这是一个值得继续深入研究的问题。
3个拟合结果显示在下图中。
使用scipy.optimize.optimize.curve_fit拟合
scipy提供的拟合，貌似需要先确定带参数的曲线方程，然后由scipy求解方程，返回曲线参数。我们还是以上面的一组数据为例使用scipy拟合曲线。
可以看出，曲线近似正弦函数。构建函数y=a*sin(x*pi/6+b)+c，使用scipy的optimize.curve_fit函数求出a、b、c的值：
求解非线性方程（组）
在数学建模中，需要对一些稀奇古怪的方程（组）求解，Matlab自然是首选，但Matlab不是免费的，scipy则为我们提供了免费的午餐！scipy.optimize库中的fsolve函数可以用来对非线性方程（组）进行求解。它的基本调用形式如下：
fsolve(func, x0)
func(x)是计算方程组误差的函数，它的参数x是一个矢量，表示方程组的各个未知数的一组可能解，func返回将x代入方程组之后得到的误差；x0为未知数矢量的初始值。
我们先来求解一个简单的方程：sin(x)−cos(x)=0.2
>>> from scipy.optimize import fsolve
>>> import numpy as np
>>> def f(A):
x = float(A[0])
return [np.sin(x) - np.cos(x) - 0.2]
>>> result = fsolve(f, [1])
array([ 0.92729522])
>>> print result
[0.92729522]
>>> print f(result)
[2.7977428707082197e-09]
哈哈，易如反掌！再来一个方程组：
4x2−2sin(yz)=0
5y+3=0
yz−1.5=0
图像处理
在scipy.misc模块中，有一个函数可以载入Lena图像——这副图像是被用作图像处理的经典示范图像。我只是简单展示一下在该图像上的几个操作。
载入Lena图像，并显示灰度图像
应用中值滤波扫描信号的每一个数据点，并替换为相邻数据点的中值
旋转图像
应用Prewitt滤波器（基于图像强度的梯度计算）
后记
这篇博文自2016年9月初动笔，断断续续写了5个多月。延宕这么久，除了自身懒惰的原因外，主要是因为MSN这个主题涉及的内容太过繁杂，又极其晦涩，无论怎么努力，总怕挂一漏万、贻笑大方。
现在好了，终于写完了。倘若哪位看官发现了谬误，请自行修改，顺便通知我一声；若因此文受益而想约饭、约酒，请发邮件至：
xufive@gmail.com
【关于作者】
许向武：山东远思信息科技有限公司CEO，网名牧码人，齐国土著，太公之后。少小离家，独闯江湖，后归隐于华不注山。素以敲击键盘为业，偶尔游戏于各网络对局室，擅长送财送分，深为众棋友所喜闻乐见。
【关于投稿】
如果大家有原创好文投稿，请直接给公号发送留言。
① 留言格式：
【投稿】+《 文章标题》+ 文章链接
② 示例：
【投稿】《不要自称是程序员，我十多年的 IT 职场总结》：http://blog.jobbole.com/94148/
③ 最后请附上您的个人简介哈~
看完本文有收获？请转发分享给更多人
关注「Python开发者」，提升Python技能