IPython 交互式计算和可视化秘籍第二版（五）
原文：
annas-archive.org/md5/62516af4e05f6a3b0523d8aa07fd5acb译者：飞龙
第十三章：随机动力学系统
在本章中，我们将讨论以下主题：
-
模拟离散时间马尔科夫链
-
模拟泊松过程
-
模拟布朗运动
-
模拟随机微分方程
介绍
随机动力学系统是受噪声影响的动力学系统。噪声带来的随机性考虑了现实世界现象中观察到的变化性。例如，股价的演变通常表现为长期行为，并伴有较快、幅度较小的振荡，反映了日常或小时的波动。
随机系统在数据科学中的应用包括统计推断方法（如马尔科夫链蒙特卡洛）和用于时间序列或地理空间数据的随机建模。
随机离散时间系统包括离散时间马尔科夫链。马尔科夫性质意味着系统在时间n+1时刻的状态仅依赖于它在时间n时刻的状态。随机元胞自动机是元胞自动机的随机扩展，是特殊的马尔科夫链。
至于连续时间系统，带噪声的常微分方程会得到随机微分方程（SDEs）。带噪声的偏微分方程会得到随机偏微分方程（SPDEs）。
点过程是另一种随机过程。这些过程建模了随时间（例如排队中顾客的到达或神经系统中的动作电位）或空间（例如森林中树木的位置、区域中的城市或天空中的星星）随机发生的瞬时事件。
从数学上讲，随机动力学系统的理论基于概率论和测度论。连续时间随机系统的研究建立在随机微积分的基础上，随机微积分是对微积分（包括导数和积分）的扩展，适用于随机过程。
在本章中，我们将看到如何使用 Python 模拟不同种类的随机系统。
参考资料
这里有一些相关的参考资料：
-
随机动力学系统概述，见于 www.scholarpedia.org/article/Sto…
-
维基百科上的马尔科夫性质，见于
en.wikipedia.org/wiki/Markov_property
模拟离散时间马尔科夫链
离散时间马尔科夫链是随机过程，它在状态空间中从一个状态转换到另一个状态。每个时间步都会发生状态转换。马尔科夫链的特点是没有记忆，即从当前状态到下一个状态的转换概率仅依赖于当前状态，而不依赖于之前的状态。这些模型在科学和工程应用中得到了广泛的使用。
连续时间马尔可夫过程也存在，我们将在本章稍后讨论特定的实例。
马尔可夫链在数学上相对容易研究，并且可以通过数值方法进行模拟。在这个案例中，我们将模拟一个简单的马尔可夫链，模拟种群的演化。
如何执行...
-
让我们导入 NumPy 和 matplotlib：
In [1]: import numpy as np import matplotlib.pyplot as plt %matplotlib inline -
我们考虑一个最大只能包含N=100个个体的种群，并定义出生和死亡率：
In [2]: N = 100 # maximum population size a = 0.5/N # birth rate b = 0.5/N # death rate -
我们在有限空间*{0, 1, ..., N}上模拟一个马尔可夫链。每个状态代表一个种群大小。
x向量将在每个时间步包含种群大小。我们将初始状态设置为x[0]=25*（即初始化时种群中有 25 个个体）：In [3]: nsteps = 1000 x = np.zeros(nsteps) x[0] = 25 -
现在我们模拟我们的链。在每个时间步t，有一个以ax[t]的概率出生，并且独立地，有一个以bx[t]的概率死亡。这些概率与当时的种群大小成正比。如果种群大小达到0或N，演化停止：
In [4]: for t in range(nsteps - 1): if 0 < x[t] < N-1: # Is there a birth? birth = np.random.rand() <= a*x[t] # Is there a death? death = np.random.rand() <= b*x[t] # We update the population size. x[t+1] = x[t] + 1*birth - 1*death # The evolution stops if we reach 0 or N. else: x[t+1] = x[t] -
让我们看看种群大小的演变：
In [5]: plt.plot(x)我们看到，在每个时间步，种群大小可以保持稳定，增加，或者减少 1。
-
现在，我们将模拟多个独立的马尔可夫链试验。我们可以用循环运行之前的仿真，但那样会非常慢（两个嵌套的
for循环）。相反，我们通过一次性考虑所有独立的试验来向量化仿真。这里只有一个关于时间的循环。在每个时间步，我们用向量化操作同时更新所有试验的结果。x向量现在包含所有试验在特定时间的种群大小。在初始化时，种群大小被设置为介于0和N之间的随机数：In [6]: ntrials = 100 x = np.random.randint(size=ntrials, low=0, high=N) -
我们定义一个执行仿真的函数。在每个时间步，我们通过生成随机向量来找到经历出生和死亡的试验，并通过向量操作更新种群大小：
In [7]: def simulate(x, nsteps): """Run the simulation.""" for _ in range(nsteps - 1): # Which trials to update? upd = (0 < x) & (x < N-1) # In which trials do births occur? birth = 1*(np.random.rand(ntrials) <= a*x) # In which trials do deaths occur? death = 1*(np.random.rand(ntrials) <= b*x) # We update the population size for all # trials. x[upd] += birth[upd] - death[upd] -
现在，让我们来看一下在不同时间点的种群大小的直方图。这些直方图代表了马尔可夫链的概率分布，通过独立试验（蒙特卡洛方法）进行估计：
In [8]: bins = np.linspace(0, N, 25) In [9]: nsteps_list = [10, 1000, 10000] for i, nsteps in enumerate(nsteps_list): plt.subplot(1, len(nsteps_list), i + 1) simulate(x, nsteps) plt.hist(x, bins=bins) plt.xlabel("Population size") if i == 0: plt.ylabel("Histogram") plt.title("{0:d} time steps".format(nsteps))而且，最初，种群大小看起来在0和N之间均匀分布，但经过足够长的时间后，它们似乎会收敛到0或N。这是因为0和N是吸收状态；一旦到达这些状态，链就无法离开这些状态。而且，这些状态可以从任何其他状态到达。
它是如何工作的...
从数学上讲，一个离散时间马尔可夫链在空间E上是一个随机变量序列X[1], X[2], ...，它满足马尔可夫性质：
一个（平稳的）马尔可夫链的特点是转移概率 P(X[j] | X[i])。这些值构成一个矩阵，称为转移矩阵。这个矩阵是一个有向图的邻接矩阵，称为状态图。每个节点是一个状态，如果链条在这些节点之间有非零转移概率，则节点 i 会连接到节点 j。
还有更多内容...
在 Python 中模拟单个马尔可夫链效率并不高，因为我们需要使用 for 循环。然而，通过向量化和并行化，可以高效地模拟许多独立的链条，这些链条遵循相同的过程（所有任务都是独立的，因此该问题是令人尴尬地并行的）。当我们关注链条的统计性质时，这种方法非常有用（例如蒙特卡罗方法的例子）。
关于马尔可夫链的文献非常广泛。许多理论结果可以通过线性代数和概率论来建立。你可以在维基百科上找到相关参考文献和教科书。
离散时间马尔可夫链有很多推广。马尔可夫链可以定义在无限状态空间上，或者在连续时间上定义。此外，马尔可夫性质在广泛的随机过程类别中非常重要。
以下是一些参考文献：
-
维基百科上的马尔可夫链，网址为
en.wikipedia.org/wiki/Markov_chain -
吸收马尔可夫链，维基百科上有相关介绍，网址为
en.wikipedia.org/wiki/Absorbing_Markov_chain -
蒙特卡罗方法，维基百科上有相关介绍，网址为
en.wikipedia.org/wiki/Monte_Carlo_method
另见
- 模拟布朗运动 示例
模拟泊松过程
泊松过程是一种特定类型的点过程，它是一种随机模型，表示瞬时事件的随机发生。大致来说，泊松过程是最不结构化或最随机的点过程。
泊松过程是一个特定的连续时间马尔可夫过程。
点过程，特别是泊松过程，可以用来模拟随机瞬时事件，比如客户在队列或服务器中的到达、电话呼叫、放射性衰变、神经细胞的动作电位等许多现象。
在这个示例中，我们将展示模拟齐次平稳泊松过程的不同方法。
如何实现...
-
让我们导入 NumPy 和 matplotlib：
In [1]: import numpy as np import matplotlib.pyplot as plt %matplotlib inline -
让我们指定
rate值，即每秒钟发生事件的平均次数：In [2]: rate = 20\. # average number of events per second -
首先，我们将使用 1 毫秒的小时间单元模拟这个过程：
In [3]: dt = .001 # time step n = int(1./dt) # number of time steps -
在每个时间区间，事件发生的概率大约是 rate * dt，如果 dt 足够小。此外，由于泊松过程没有记忆性，事件的发生在不同的区间之间是独立的。因此，我们可以在向量化的方式中采样伯努利随机变量（分别表示实验的成功或失败，即 1 或 0），以此来模拟我们的过程：
In [4]: x = np.zeros(n) x[np.random.rand(n) <= rate*dt] = 1x向量包含了所有时间区间上的零和一，1表示事件的发生：In [5]: x[:5] Out[5]: array([ 0., 1., 0., 0., 0\. ]) -
我们来展示模拟过程。我们为每个事件绘制一条垂直线：
In [6]: plt.vlines(np.nonzero(x)[0], 0, 1) plt.xticks([]); plt.yticks([]) -
表示同一对象的另一种方式是考虑相关的 计数过程 N(t)，它表示在时间 t 之前发生的事件数。在这里，我们可以使用
cumsum()函数来展示这个过程：In [7]: plt.plot(np.linspace(0., 1., n), np.cumsum(x)) plt.xlabel("Time") plt.ylabel("Counting process") -
模拟均匀泊松过程的另一种（更高效的）方式是利用两个连续事件之间的时间间隔遵循指数分布这一性质。此外，这些间隔是独立的。因此，我们可以以向量化的方式对它们进行采样。最终，通过累加这些间隔，我们得到了我们的过程：
In [8]: y = np.cumsum(np.random.exponential( 1./rate, size=int(rate)))y向量包含了我们的泊松过程的另一个实现，但数据结构不同。向量的每个分量都是一个事件发生的时间：In [9]: y[:5] Out[9]: array([ 0.006, 0.111, 0.244, 0.367, 0.365]) -
最后，我们来展示模拟过程：
In [10]: plt.vlines(y, 0, 1) plt.xticks([]); plt.yticks([])
它是如何工作的...
对于一个速率为 的泊松过程，长度为
的时间窗口内的事件数服从泊松分布：
当 很小时，我们可以证明，按照一阶近似，这个概率大约是
。
此外，保持时间（两个连续事件之间的延迟）是独立的，并且遵循指数分布。泊松过程满足其他有用的性质，如独立和稳定增量。这一性质证明了本方法中使用的第一次模拟方法。
还有更多...
在这个示例中，我们只考虑了均匀时间相关的泊松过程。其他类型的泊松过程包括不均匀（或非均匀）过程，其特点是速率随时间变化，以及多维空间泊松过程。
这里有更多的参考资料：
-
维基百科上的泊松过程，链接：
en.wikipedia.org/wiki/Poisson_process -
维基百科上的点过程，链接：
en.wikipedia.org/wiki/Point_process -
维基百科上的连续时间过程，链接：
en.wikipedia.org/wiki/Continuous-time_process -
维基百科上的更新理论，可在
en.wikipedia.org/wiki/Renewal_theory中查看 -
维基百科上的空间泊松过程，可在
en.wikipedia.org/wiki/Spatial_Poisson_process中查看
另请参见
- 模拟离散时间马尔可夫链方法
模拟布朗运动
布朗运动（或维纳过程）是数学、物理学以及许多其他科学和工程学科中的基本对象。该模型描述了粒子在流体中由于与流体中快速分子发生随机碰撞而产生的运动（扩散）。更一般地，布朗运动模型描述了一种连续时间随机游走，其中粒子通过在所有方向上做独立的随机步骤在空间中演化。
从数学角度来看，布朗运动是一个特殊的马尔可夫连续随机过程。布朗运动是随机微积分和随机过程理论等数学领域的核心，但在量化金融、生态学和神经科学等应用领域也占有重要地位。
在这个方法中，我们将展示如何在二维空间中模拟和绘制布朗运动。
如何实现...
-
让我们导入 NumPy 和 matplotlib：
In [1]: import numpy as np import matplotlib.pyplot as plt %matplotlib inline -
我们用 5000 个时间步骤来模拟布朗运动：
In [2]: n = 5000 -
我们模拟两个独立的一维布朗过程来形成一个二维布朗过程。（离散）布朗运动在每个时间步骤上做独立的高斯跳跃。因此，我们只需要计算独立正态随机变量的累积和（每个时间步骤一个）：
In [3]: x = np.cumsum(np.random.randn(n)) y = np.cumsum(np.random.randn(n)) -
现在，为了展示布朗运动，我们可以直接使用
plot(x, y)。然而，这样的结果会是单色的，稍显乏味。我们希望使用渐变色来展示运动随时间的进展（色调是时间的函数）。matplotlib 迫使我们使用基于scatter()的小技巧。这个函数允许我们为每个点分配不同的颜色，代价是丢失了点与点之间的线段。为了解决这个问题，我们线性插值这个过程，以便呈现出连续线条的错觉：In [4]: k = 10 # We add 10 intermediary points between two # successive points. # We interpolate x and y. x2 = np.interp(np.arange(n*k), np.arange(n)*k, x) y2 = np.interp(np.arange(n*k), np.arange(n)*k, y) In [5]: # Now, we draw our points with a gradient of # colors. plt.scatter(x2, y2, c=range(n*k), linewidths=0, marker='o', s=3, cmap=plt.cm.jet) plt.axis('equal') plt.xticks([]); plt.yticks([])
如何实现...
布朗运动*W(t)*具有几个重要的性质。首先，它几乎肯定会产生连续的轨迹。其次，它的增量在不重叠的区间上是独立的。第三，这些增量是高斯随机变量。更准确地说：
特别地，W(t)的密度是一个方差为t的正态分布。
此外，布朗运动以及一般的随机过程与偏微分方程有着深刻的联系。在这里，W(t) 的密度是 热方程 的解，这是一种特殊的扩散方程。更一般地说，Fokker-Planck 方程 是由随机微分方程的解的密度所满足的偏微分方程。
还有更多...
布朗运动是具有无限小步长的随机游走的极限。在这里，我们利用这一特性来模拟该过程。
这里有一些参考文献：
-
在
en.wikipedia.org/wiki/Brownian_motion上描述了布朗运动（物理现象） -
在
en.wikipedia.org/wiki/Wiener_process上解释了维纳过程（数学对象） -
布朗运动是 Lévy 过程的一种特例；参见
en.wikipedia.org/wiki/L%C3%A9vy_process -
Fokker-Planck 方程将随机过程与偏微分方程联系起来；参见
en.wikipedia.org/wiki/Fokker%E2%80%93Planck_equation
另见
- 模拟随机微分方程的步骤
模拟随机微分方程
随机微分方程（SDEs）模拟受噪声影响的动态系统。它们广泛应用于物理学、生物学、金融学以及其他学科。
在本步骤中，我们模拟一个 Ornstein-Uhlenbeck 过程，它是 Langevin 方程的解。该模型描述了在摩擦作用下粒子在流体中的随机演化。粒子的运动是由于与流体分子发生碰撞（扩散）。与布朗运动的区别在于存在摩擦。
Ornstein-Uhlenbeck 过程是平稳的、高斯的和马尔可夫的，这使得它成为表示平稳随机噪声的良好候选模型。
我们将通过一种叫做 Euler-Maruyama 方法 的数值方法来模拟这个过程。它是对常微分方程（ODE）中欧拉方法的一个简单推广，适用于随机微分方程（SDE）。
如何做...
-
让我们导入 NumPy 和 matplotlib：
In [1]: import numpy as np import matplotlib.pyplot as plt %matplotlib inline -
我们为模型定义了一些参数：
In [2]: sigma = 1\. # Standard deviation. mu = 10.0 # Mean. tau = 0.05 # Time constant. -
让我们定义一些模拟参数：
In [3]: dt = 0.001 # Time step. T = 1.0 # Total time. n = int(T/dt) # Number of time steps. t = np.linspace(0., T, n) # Vector of times. -
我们还定义了重正规化变量（以避免在每个时间步长重新计算这些常数）：
In [4]: sigma_bis = sigma * np.sqrt(2\. / tau) sqrtdt = np.sqrt(dt) -
我们创建一个向量，在模拟过程中包含所有连续的过程值：
In [5]: x = np.zeros(n) -
现在，让我们使用 Euler-Maruyama 方法来模拟该过程。它实际上是常规欧拉法在常微分方程（ODE）中的一种推广，但增加了一个额外的随机项（它仅仅是一个缩放的标准正态随机变量）。我们将在 如何工作... 部分给出该过程的方程和该方法的详细信息：
In [6]: for i in range(n-1): x[i+1] = x[i] + dt*(-(x[i]-mu)/tau) + \ sigma_bis * sqrtdt * np.random.randn() -
让我们展示过程的演变：
In [7]: plt.plot(t, x) -
现在，我们将查看该过程分布随时间演化的情况。为此，我们将以向量化方式模拟该过程的多个独立实现。我们定义一个向量
X，它将包含在给定时间点所有实现的过程（也就是说，我们不会在所有时间点都将所有实现保留在内存中）。该向量将在每个时间步覆盖。我们将在多个时间点展示估计的分布（直方图）：In [8]: ntrials = 10000 X = np.zeros(ntrials) In [9]: # We create bins for the histograms. bins = np.linspace(-2., 14., 50); for i in range(n): # We update the process independently for all # trials. X += dt*(-(X-mu)/tau) + \ sigma_bis*sqrtdt*np.random.randn(ntrials) # We display the histogram for a few points in # time. if i in (5, 50, 900): hist, _ = np.histogram(X, bins=bins) plt.plot((bins[1:]+bins[:-1])/2, hist, label="t={0:.2f}".format(i*dt)) plt.legend()该过程的分布趋近于均值为
和标准差为
的高斯分布。如果初始分布也是具有适当参数的高斯分布，那么该过程将是平稳的。
它是如何工作的...
我们在这个实例中使用的 Langevin 方程是以下随机微分方程：
在这里，x(t) 是我们的随机过程，dx 是无穷小增量， 是均值，
是标准差，
是时间常数。此外，W 是布朗运动（或维纳过程），它是我们随机微分方程的基础。
右边的第一个项是确定性项（在 dt 中），而第二个项是随机项。如果没有最后一个项，该方程将是一个常规的确定性常微分方程。
布朗运动的无穷小步长是一个高斯随机变量。具体来说，布朗运动的导数（以某种意义上）是 白噪声，即一系列独立的高斯随机变量。
欧拉-马鲁亚马方法涉及时间离散化，并在每个时间步向过程添加无穷小步长。该方法包含一个确定性项（如常规欧拉方法在常微分方程中的应用）和一个随机项（随机高斯变量）。具体来说，对于一个方程：
数值方案是（其中 t=n * dt）：
在这里， 是一个方差为 1 的随机高斯变量（在每个时间步独立）。归一化因子
来源于布朗运动的无穷小步长，其标准差为
。
还有更多...
随机微分方程的数学理论包括随机微积分、伊藤微积分、马尔可夫过程等主题。尽管这些理论相当复杂，但正如我们在这个实例中看到的那样，数值模拟随机过程相对简单。
欧拉-马鲁亚马方法的误差是 dt 量级的。Milstein 方法是一种更精确的数值方案，误差量级为 dt。
这里有一些相关参考资料：
-
维基百科上的随机微分方程
-
维基百科上的朗之万方程
-
奥恩斯坦-乌伦贝克过程的描述
-
欧拉-马鲁亚马方法，解释于
en.wikipedia.org/wiki/Euler%E2%80%93Maruyama_method -
维基百科上的米尔斯坦方法
另请参见
- 模拟布朗运动的方案
第十四章：图、几何学和地理信息系统
在本章中，我们将涵盖以下主题：
-
使用 NetworkX 操作和可视化图
-
使用 NetworkX 分析社交网络
-
使用拓扑排序解决有向无环图中的依赖关系
-
计算图像中的连通分量
-
计算一组点的沃罗诺伊图
-
使用 Shapely 和 basemap 操作地理空间数据
-
为路网创建路径规划器
介绍
在本章中，我们将涵盖 Python 在图论、几何学和地理学方面的能力。
图是描述物体之间关系的数学对象。它们在科学和工程中无处不在，因为它们可以表示许多现实世界的关系：社交网络中的朋友、分子中的原子、网站链接、神经网络中的细胞、图像中的相邻像素等。图还是计算机科学中的经典数据结构。最后，许多领域特定的问题可以重新表达为图论问题，然后使用著名的算法来解决。
我们还将看到一些与几何和地理信息系统（GIS）相关的内容，GIS 指的是对任何类型的空间、地理或地形数据的处理和分析。
在本介绍中，我们将简要概述这些主题。
图
从数学角度来看，图 G = (V, E) 由一组顶点或节点 V 和一组边 E（V 的二元子集）定义。若 (v, v') 是一条边（E 的元素），则称两个节点 v 和 v' 是连接的。
-
如果边是无序的（意味着 (v,v') = (v',v)），则图被称为无向图
-
如果边是有序的（意味着 (v,v') ≠ (v',v)），则图被称为有向图
在无向图中，边由连接两个节点的线段表示。在有向图中，它由箭头表示。
无向图和有向图
图可以通过不同的数据结构表示，特别是邻接表（每个顶点的邻接顶点列表）或邻接矩阵（顶点之间连接的矩阵）。
图论中的问题
这里有几个经典的图论问题的示例：
-
图遍历：如何遍历图，详细讨论见
en.wikipedia.org/wiki/Graph_traversal -
图着色：如何给图中的节点着色，使得相邻的两个顶点不共享相同的颜色，详细讨论见
en.wikipedia.org/wiki/Graph_coloring -
连通分量：如何在图中找到连通分量，详细解释见
en.wikipedia.org/wiki/Connected_component_%28graph_theory%29 -
最短路径：在给定图中，从一个节点到另一个节点的最短路径是什么？讨论见
en.wikipedia.org/wiki/Shortest_path_problem -
哈密尔顿路径：一个图是否包含哈密尔顿路径，访问每个顶点恰好一次？详细解释见
en.wikipedia.org/wiki/Hamiltonian_path -
欧拉路径：一个图是否包含欧拉路径，访问每条 边 恰好一次？讨论见
en.wikipedia.org/wiki/Eulerian_path -
旅行商问题：访问每个节点恰好一次（哈密尔顿路径）的最短路径是什么？详细解释见
en.wikipedia.org/wiki/Traveling_salesman_problem
随机图
随机图 是一种特殊的图，通过概率规则定义。它们有助于理解大规模现实世界图的结构，如社交图。
特别地，小世界网络具有稀疏的连接，但大多数节点可以在少数几步内从其他任何节点到达。这一特性源于少数几个具有大量连接的 枢纽 的存在。
Python 中的图
尽管可以使用原生的 Python 结构操作图，但使用专门的库来实现特定数据结构和操作例程会更方便。在本章中，我们将使用 NetworkX，一个纯 Python 库。其他 Python 库包括 python-graph 和 graph-tool（主要用 C++ 编写）。
NetworkX 实现了一个灵活的图数据结构，并包含许多算法。NetworkX 还允许我们使用 matplotlib 轻松绘制图。
Python 中的几何学
Shapely 是一个 Python 库，用于处理二维几何图形，如点、线和多边形。它在地理信息系统中尤其有用。
将 Shapely 与 matplotlib 结合并非易事。幸运的是，descartes 包使这个任务变得更加简单。
Python 中的地理信息系统（GIS）
有几个 Python 模块用于处理地理数据和绘制地图。
在本章中，我们将使用 matplotlib 的 basemap、Shapely、descartes 和 Fiona 来处理 GIS 文件。
ESRI Shapefile 是一种流行的地理空间矢量数据格式。它可以被 basemap、NetworkX 和 Fiona 读取。
我们还将使用 OpenStreetMap 服务，这是一个免费的开源协作服务，提供全球地图。
本章未涉及的其他 Python 中的 GIS/地图系统包括 GeoPandas、Kartograph、Vincent 和 cartopy。
参考资料
以下是关于图的一些参考资料：
-
维基百科中的图论，访问
en.wikipedia.org/wiki/Graph_theory -
描述图的数据结构，请参考
en.wikipedia.org/wiki/Graph_(abstract_data_type) -
可在维基百科上查看随机图的页面，网址为
en.wikipedia.org/wiki/Random_graph -
可在维基百科上查看小世界图的页面，网址为
en.wikipedia.org/wiki/Small-world_network -
NetworkX 软件包的网址为
networkx.github.io -
python-graph 软件包的网址为
code.google.com/p/python-graph/ -
可以在
graph-tool.skewed.de获取 graph-tool 软件包。
以下是关于 Python 中几何和地图的一些参考资料：
-
Basemap 的网址为
matplotlib.org/basemap/ -
Shapely 的网址为
toblerity.org/shapely/project.html -
Fiona 的网址为
toblerity.org/fiona/ -
descartes 的网址为
pypi.python.org/pypi/descartes -
Shapefile 的网址为
en.wikipedia.org/wiki/Shapefile -
OpenStreetMap 的网址为www.openstreetmap.org
-
Folium 的网址为
github.com/wrobstory/folium -
GeoPandas 的网址为
geopandas.org -
Kartograph 的网址为
kartograph.org -
Cartopy 的网址为
scitools.org.uk/cartopy/ -
Vincent 的网址为
github.com/wrobstory/vincent
使用 NetworkX 操作和可视化图形
在这个示例中，我们将展示如何使用 NetworkX 创建、操作和可视化图形。
准备工作
您可以在官方文档中找到 NetworkX 的安装说明，网址为networkx.github.io/documentation/latest/install.html。
使用 Anaconda，您可以在终端中输入conda install networkx。或者，您也可以输入pip install networkx。在 Windows 上，您还可以使用 Chris Gohlke 的安装程序，网址为www.lfd.uci.edu/~gohlke/pyt…。
如何做…
-
让我们导入 NumPy、NetworkX 和 matplotlib：
In [1]: import numpy as np import networkx as nx import matplotlib.pyplot as plt %matplotlib inline -
创建图形的方法有很多种。在这里，我们创建了一个边的列表（节点索引的对）：
In [2]: n = 10 # Number of nodes in the graph. # Each node is connected to the two next nodes, # in a circular fashion. adj = [(i, (i+1)%n) for i in range(n)] adj += [(i, (i+2)%n) for i in range(n)] -
我们使用边的列表实例化一个
Graph对象：In [3]: g = nx.Graph(adj) -
让我们检查图的节点和边的列表，以及其邻接矩阵：
In [4]: print(g.nodes()) [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] In [5]: print(g.edges()) [(0, 8), (0, 1), (0, 2), ..., (7, 9), (8, 9)] In [6]: print(nx.adjacency_matrix(g)) [[ 0\. 1\. 1\. 0\. 0\. 0\. 0\. 0\. 1\. 1.] [ 1\. 0\. 1\. 1\. 0\. 0\. 0\. 0\. 0\. 1.] ... [ 1\. 0\. 0\. 0\. 0\. 0\. 1\. 1\. 0\. 1.] [ 1\. 1\. 0\. 0\. 0\. 0\. 0\. 1\. 1\. 0.]] -
让我们显示这个图。NetworkX 自带各种绘图函数。我们可以明确指定节点的位置，也可以使用算法自动计算一个有趣的布局。在这里，我们使用了
draw_circular()函数，它会将节点线性地放置在一个圆上：In [7]: nx.draw_circular(g) -
图可以很容易地修改。在这里，我们添加一个连接到所有现有节点的新节点。我们还为此节点指定了一个
color属性。在 NetworkX 中，每个节点和边都带有一个方便的 Python 字典，包含任意属性。In [8]: g.add_node(n, color='#fcff00') # We add an edge from every existing # node to the new node. for i in range(n): g.add_edge(i, n) -
现在，让我们再次绘制修改后的图。这次，我们明确指定节点的位置和颜色：
In [9]: # We define custom node positions on a circle # except the last node which is at the center. t = np.linspace(0., 2*np.pi, n) pos = np.zeros((n+1, 2)) pos[:n,0] = np.cos(t) pos[:n,1] = np.sin(t) # A node's color is specified by its 'color' # attribute, or a default color if this attribute # doesn't exist. color = [g.node[i].get('color', '#88b0f3') for i in range(n+1)] # We now draw the graph with matplotlib. nx.draw_networkx(g, pos=pos, node_color=color) plt.axis('off') -
让我们也使用自动布局算法：
In [10]: nx.draw_spectral(g, node_color=color) plt.axis('off')
还有更多…
在 NetworkX 中，节点不一定是整数。它们可以是数字、字符串、元组和任何可哈希的 Python 类的实例。
此外，每个节点和边都带有可选属性（形成一个字典）。
在 NetworkX 中实现了一些布局算法。draw_spectral()函数使用图的拉普拉斯矩阵的特征向量。
draw_spring()函数实现了Fruchterman-Reingold 力导向算法。节点被视为受边缘相关力的质点。力导向图绘制算法通过最小化系统的能量来找到平衡配置。这将导致一个美观的布局，尽可能少地交叉边。
这里有一些参考资料：
-
维基百科上的拉普拉斯矩阵，网址为
en.wikipedia.org/wiki/Laplacian_matrix -
描述在
en.wikipedia.org/wiki/Force-directed_graph_drawing的力导向图绘制
另请参阅
- 使用 NetworkX 分析社交网络配方
使用 NetworkX 分析社交网络
在这个配方中，我们将展示如何在 Python 中分析社交数据。社交数据是由人们在社交网络上的活动生成的，如 Facebook、Twitter、Google+、GitHub 等。
在这个配方中，我们将使用 NetworkX 分析和可视化 Twitter 用户的社交网络。
准备工作
首先，您需要安装Twitter Python 包。您可以使用pip install twitter进行安装。您可以在pypi.python.org/pypi/twitter找到更多信息。
然后，您需要获取认证代码以访问您的 Twitter 数据。该过程是免费的。除了 Twitter 账号外，您还需要在 Twitter 开发者网站上创建一个应用程序，网址为dev.twitter.com/apps。然后，您将能够检索到此配方所需的OAuth 认证代码。
您需要在当前文件夹中创建一个twitter.txt文本文件，其中包含四个私有认证密钥。每行必须有一个密钥，顺序如下：
-
API 密钥
-
API 秘钥
-
访问令牌
-
访问令牌密钥
请注意，对 Twitter API 的访问是有限制的。大多数方法在给定时间窗口内只能调用几次。除非您研究小网络或查看大网络的小部分，否则您将需要节流您的请求。在这个示例中，我们只考虑网络的一小部分，因此不应该达到 API 限制。否则，您将需要等待几分钟，直到下一个时间窗口开始。API 限制可在dev.twitter.com/docs/rate-limiting/1.1/limits查看。
如何做…
-
让我们导入一些包：
In [1]: import math import json import twitter import numpy as np import pandas as pd import networkx as nx import matplotlib.pyplot as plt %matplotlib inline from IPython.display import Image -
我们从我们的
twitter.txt文件中获取秘密的消费者和 OAuth 密钥：In [2]: (CONSUMER_KEY, CONSUMER_SECRET, OAUTH_TOKEN, OAUTH_TOKEN_SECRET) = open( 'twitter.txt', 'r').read().splitlines() -
现在我们创建一个
Twitter类的实例，这将使我们能够访问 Twitter API：In [3]: auth = twitter.oauth.OAuth(OAUTH_TOKEN, OAUTH_TOKEN_SECRET, CONSUMER_KEY, CONSUMER_SECRET) tw = twitter.Twitter(auth=auth) -
在这个示例中，我们使用 Twitter API 的 1.1 版本。
twitter库在Twitter实例的属性之间定义了 REST API 和直接映射。在这里，我们执行account/verify_credentialsREST 请求来获取认证用户的标识符（这里是我，或者如果您自己执行这个笔记本，则是您！）：In [4]: me = tw.account.verify_credentials() In [5]: myid = me['id'] -
让我们定义一个简单的函数，返回给定用户（默认为认证用户）的所有关注者的标识符：
In [6]: def get_followers_ids(uid=None): # Retrieve the list of followers' ids of the # specified user. return tw.followers.ids(user_id=uid)['ids'] In [7]: # We get the list of my followers. my_followers_ids = get_followers_ids() -
现在，我们定义一个函数，用于检索 Twitter 用户的完整资料。由于
users/lookup批量请求每次限制为 100 个用户，并且在时间窗口内只允许少量调用，我们只查看所有关注者的一个子集：In [8]: def get_users_info(users_ids, max=500): n = min(max, len(users_ids)) # Get information about those users, # using batch requests. users = [tw.users.lookup( user_id=users_ids[100*i:100*(i+1)]) for i in range(int(math.ceil(n/100.)))] # We flatten this list of lists. users = [item for sublist in users for item in sublist] return {user['id']: user for user in users} In [9]: users_info = get_users_info(my_followers_ids) In [10]: # Let's save this dictionary on the disk. with open('my_followers.json', 'w') as f: json.dump(users_info, f, indent=1) -
我们还开始用关注者定义图形，使用邻接表（技术上来说，是一个列表的字典）。这被称为自我图。这个图表示我们的关注者之间的所有关注连接：
In [11]: adjacency = {myid: my_followers_ids} -
现在，我们将查看与 Python 相关的自我图部分。具体来说，我们将考虑包含“Python”描述的前 10 位最受关注用户的关注者：
In [12]: my_followers_python = \[user for user in users_info.values() if 'python' in user['description'].lower()] In [13]: my_followers_python_best = \ sorted(my_followers_python, key=lambda u: u['followers_count'])[::-1][:10]检索给定用户的关注者的请求受到速率限制。让我们来看看我们还剩下多少次调用：
In [14]: tw.application.rate_limit_status( resources='followers') \ ['resources']['followers']['/followers/ids'] Out[14]: {u'limit': 15, u'remaining': 0, u'reset': 1388865551} In [15]: for user in my_followers_python_best: # The call to get_followers_ids is # rate-limited. adjacency[user['id']] = list(set( get_followers_ids(user['id'])). \ intersection(my_followers_ids)) -
现在我们的图被定义为字典中的邻接表，我们将在 NetworkX 中加载它：
In [16]: g = nx.Graph(adjacency) In [17]: # We only restrict the graph to the users # for which we were able to retrieve the profile. g = g.subgraph(users_info.keys()) In [18]: # We also save this graph on disk. with open('my_graph.json', 'w') as f: json.dump(nx.to_dict_of_lists(g), f, indent=1) In [19]: # We remove isolated nodes for simplicity. g.remove_nodes_from([k for k, d in g.degree().items() if d <= 1]) In [20]: # Since I am connected to all nodes, # by definition, we can remove me for simplicity. g.remove_nodes_from([myid]) -
让我们看一下图的统计数据：
In [21]: len(g.nodes()), len(g.edges()) Out[21]: (197, 1037) -
现在我们将绘制这个图。我们将根据每个用户的关注者数量和推文数量使用不同的大小和颜色来绘制节点。最受关注的用户将更大。最活跃的用户将更红。
In [22]: # Update the dictionary. deg = g.degree() for user in users_info.values(): fc = user['followers_count'] sc = user['statuses_count'] # Is this user a Pythonista? user['python'] = 'python' in \ user['description'].lower() # We compute the node size as a function of # the number of followers. user['node_size'] = math.sqrt(1 + 10 * fc) # The color is function of its activity user['node_color'] = 10 * math.sqrt(1 + sc) # We only display the name of the most # followed users. user['label'] = user['screen_name'] \ if fc > 2000 else '' -
最后，我们使用
draw()函数来显示图形。我们需要将节点的大小和颜色指定为列表，标签指定为字典：In [23]: node_size = [users_info[uid]['node_size'] for uid in g.nodes()] In [24]: node_color = [users_info[uid]['node_color'] for uid in g.nodes()] In [25]: labels = {uid: users_info[uid]['label'] for uid in g.nodes()} In [26]: nx.draw(g, cmap=plt.cm.OrRd, alpha=.8, node_size=node_size, node_color=node_color, labels=labels, font_size=4, width=.1)
还有更多…
一本关于使用 Python 进行社交数据分析的重要参考资料是 Matthew A. Russel 的书籍*《Mining the Social Web》，由O'Reilly Media*出版。 代码可以在 GitHub 上的 IPython 笔记本中找到，网址为github.com/ptwobrussell/Mining-the-Social-Web-2nd-Edition。 以下网络被涵盖：Twitter、Facebook、LinkedIn、Google+、GitHub、邮箱、网站等。
参见
- 使用 NetworkX 操纵和可视化图形示例
使用拓扑排序在有向无环图中解析依赖关系
在这个示例中，我们将展示一个著名的图算法应用：拓扑排序。 让我们考虑描述项目之间依赖关系的有向图。 例如，在软件包管理器中，在安装给定软件包P之前，我们可能需要安装依赖软件包。
依赖关系集合形成一个有向图。 通过拓扑排序，软件包管理器可以解析依赖关系并找到软件包的正确安装顺序。
拓扑排序有许多其他应用。 在这里，我们将在 Debian 软件包管理器的真实数据上说明这个概念。 我们将找到 IPython 所需软件包的安装顺序。
准备工作
您需要python-apt软件包才能构建软件包依赖关系图。 该软件包可在pypi.python.org/pypi/python-apt/上找到。
我们还假设此笔记本在 Debian 系统（如 Ubuntu）上执行。 如果您没有这样的系统，可以直接从书的 GitHub 存储库github.com/ipython-books/cookbook-data下载Debian数据集。 将其解压缩到当前目录，并直接从此笔记本的第 7 步开始。
如何做…
-
我们导入
apt模块并构建软件包列表：In [1]: import json import apt cache = apt.Cache() -
graph字典将包含依赖关系图的一小部分的邻接表：In [2]: graph = {} -
我们定义一个返回软件包依赖关系列表的函数：
In [3]: def get_dependencies(package): if package not in cache: return [] pack = cache[package] ver = pack.candidate or pack.versions[0] # We flatten the list of dependencies, # and we remove the duplicates. return sorted(set([item.name for sublist in ver.dependencies for item in sublist])) -
现在我们定义一个递归函数，用于构建特定软件包的依赖关系图。 此函数更新
graph变量：In [4]: def get_dep_recursive(package): if package not in cache: return [] if package not in graph: dep = get_dependencies(package) graph[package] = dep for dep in graph[package]: if dep not in graph: graph[dep] = get_dep_recursive(dep) return graph[package] -
让我们为 IPython 构建依赖关系图：
In [5]: get_dep_recursive('ipython') -
最后，我们将邻接表保存为 JSON：
In [6]: with open('data/apt.json', 'w') as f: json.dump(graph, f, indent=1)提示
如果您没有 Debian 操作系统，请从书籍的存储库中下载Debian数据集。
-
我们导入一些包：
In [7]: import json import numpy as np import networkx as nx import matplotlib.pyplot as plt %matplotlib inline -
让我们从 JSON 文件中加载邻接表：
In [8]: with open('data/apt.json', 'r') as f: graph = json.load(f) -
现在，我们从邻接表创建一个有向图（NetworkX 中的
DiGraph）。 我们反转图以获得更自然的顺序：In [9]: g = nx.DiGraph(graph).reverse() -
当图形是有向无环图（DAG）时才存在拓扑排序。 这意味着图中没有循环，也就是说，没有循环依赖。 我们的图是 DAG 吗？ 让我们看看：
In [10]: nx.is_directed_acyclic_graph(g) Out[10]: False -
哪些软件包负责循环？ 我们可以使用
simple_cycles()函数找到它们：In [11]: set([cycle[0] for cycle in nx.simple_cycles(g)]) Out[11]: {u'coreutils', u'libc6', u'multiarch-support', u'python-minimal', u'tzdata'} -
在这里，我们可以尝试移除这些包。在实际的包管理器中，这些循环需要仔细考虑。
In [12]: g.remove_nodes_from(_) In [13]: nx.is_directed_acyclic_graph(g) Out[13]: True -
现在，图是一个 DAG。让我们先展示它：
In [14]: ug = g.to_undirected() deg = ug.degree() In [15]: # The size of the nodes depends on the number # of dependencies. nx.draw(ug, font_size=6, node_size=[20*deg[k] for k in ug.nodes()]) -
最后，我们可以执行拓扑排序，从而获得一个满足所有依赖关系的线性安装顺序：
In [16]: nx.topological_sort(g) Out[16]: [u'libexpat1', u'libdb5.1', u'debconf-2.0', ... u'python-pexpect', u'python-configobj', u'ipython']
还有更多内容…
有向无环图在许多应用中都有出现。它们可以表示因果关系、影响图、依赖关系以及其他概念。例如，像 Git 这样的分布式版本控制系统的版本历史就是用 DAG 来描述的。
拓扑排序在一般的调度任务中非常有用（项目管理和指令调度）。
这里有一些参考资料：
-
维基百科上的拓扑排序，链接：
en.wikipedia.org/wiki/Topological_sorting
计算图像中的连通分量
在本食谱中，我们将展示图论在图像处理中的应用。我们将计算图像中的连通分量。这种方法将允许我们标记图像中的连续区域，类似于绘图程序中的桶形填充工具。
找到连通分量在许多益智类视频游戏中也很有用，如扫雷、泡泡龙等。在这些游戏中，需要自动检测同一颜色的连续物品集。
如何做到这一点…
-
让我们导入这些包：
In [1]: import itertools import numpy as np import networkx as nx import matplotlib.colors as col import matplotlib.pyplot as plt %matplotlib inline -
我们创建一个10 x 10的图像，每个像素可以取三个可能的标签（或颜色）之一：
In [2]: n = 10 In [3]: img = np.random.randint(size=(n, n), low=0, high=3) -
现在，我们创建一个底层的二维网格图，编码图像的结构。每个节点是一个像素，且一个节点与其最近的邻居相连。NetworkX 定义了一个
grid_2d_graph函数来生成这个图：In [4]: g = nx.grid_2d_graph(n, n) -
让我们创建两个函数来显示图像和相应的图形：
In [5]: def show_image(img, **kwargs): plt.imshow(img, origin='lower', interpolation='none', **kwargs) plt.axis('off') In [6]: def show_graph(g, **kwargs): nx.draw(g, pos={(i, j): (j, i) for (i, j) in g.nodes()}, node_color=[img[i, j] for (i, j) in g.nodes()], linewidths=1, edge_color='w', with_labels=False, node_size=30, **kwargs) In [7]: cmap = plt.cm.Blues -
这是原始图像与底层图叠加后的效果：
In [8]: show_image(img, cmap=cmap, vmin=-1) show_graph(g, cmap=cmap, vmin=-1) -
现在，我们将找出所有包含超过三个像素的连续深蓝色区域。首先，我们考虑对应所有深蓝色像素的子图：
In [9]: g2 = g.subgraph(zip(*np.nonzero(img==2))) In [10]: show_image(img, cmap=cmap, vmin=-1) show_graph(g2, cmap=cmap, vmin=-1) -
我们看到，所请求的连续区域对应于包含超过三个节点的连通分量。我们可以使用 NetworkX 的
connected_components函数来找到这些分量：In [11]: components = [np.array(comp) for comp in nx.connected_components(g2) if len(comp)>=3] len(components) Out[11]: 3 -
最后，我们为每个分量分配一个新颜色，并显示新的图像：
In [12]: # We copy the image, and assign a new label # to each found component. img_bis = img.copy() for i, comp in enumerate(components): img_bis[comp[:,0], comp[:,1]] = i + 3 In [13]: # We create a new discrete color map extending # the previous map with new colors. colors = [cmap(.5), cmap(.75), cmap(1.), '#f4f235', '#f4a535', '#f44b35'] cmap2 = col.ListedColormap(colors, 'indexed') In [14]: show_image(img_bis, cmap=cmap2)
它是如何工作的…
我们解决的问题被称为连通分量标记。它也与洪水填充算法密切相关。
将网格图与图像关联的想法在图像处理领域非常常见。在这里，连续的颜色区域对应于连通组件的子图。一个连通组件可以定义为可达性关系的等价类。如果两个节点之间有一条路径，则它们在图中是连通的。一个等价类包含可以互相到达的节点。
最后，这里描述的简单方法仅适用于小图像上的基本任务。更高级的算法将在第十一章，图像与音频处理中讲解。
还有更多…
这里有一些参考资料：
-
Wikipedia 上的连通组件，访问
en.wikipedia.org/wiki/Connected_component_%28graph_theory%29 -
Wikipedia 上的连通组件标记算法，访问
en.wikipedia.org/wiki/Connected-component_labeling -
Wikipedia 上的 Flood-fill 算法，访问
en.wikipedia.org/wiki/Flood_fill
计算一组点的 Voronoi 图
一组种子点的Voronoi 图将空间划分为多个区域。每个区域包含所有离某个种子点比任何其他种子点更近的点。
Voronoi 图是计算几何中的基本结构。它广泛应用于计算机科学、机器人学、地理学和其他学科。例如，一组地铁站的 Voronoi 图可以告诉我们城市中任意一个点离哪个地铁站最近。
在这个示例中，我们使用 SciPy 计算巴黎地铁站集的 Voronoi 图。
准备工作
你需要 Smopy 模块来显示巴黎的 OpenStreetMap 地图。你可以使用pip install smopy安装该包。
你还需要从书籍的 GitHub 仓库下载RATP数据集，网址为github.com/ipython-books/cookbook-data，并将其解压到当前目录。数据来自 RATP 的开放数据网站（巴黎的公共交通运营商，data.ratp.fr）。
如何实现…
-
让我们导入 NumPy、pandas、matplotlib 和 SciPy：
In [1]: import numpy as np import pandas as pd import scipy.spatial as spatial import matplotlib.pyplot as plt import matplotlib.path as path import matplotlib as mpl import smopy %matplotlib inline -
让我们使用 pandas 加载数据集：
In [2]: df = pd.read_csv('data/ratp.csv', sep='#', header=None) In [3]: df[df.columns[1:]].tail(2) Out[3]: 1 2 3 4 5 11609 2.30 48.93 TIMBAUD GENNEVILLIERS tram 11610 2.23 48.91 VICTOR BASCH COLOMBES tram -
DataFrame对象包含站点的坐标、名称、城市、区域和类型。让我们选择所有地铁站：In [4]: metro = df[(df[5] == 'metro')] In [5]: metro[metro.columns[1:]].tail(3) Out[5]: 305 2.308 48.841 Volontaires PARIS-15EME metro 306 2.379 48.857 Voltaire PARIS-11EME metro 307 2.304 48.883 Wagram PARIS-17EME metro -
我们将提取巴黎各地铁站的区域编号。使用 pandas 时，我们可以通过相应列的
str属性进行向量化的字符串操作。In [6]: # We only extract the district from stations # in Paris. paris = metro[4].str.startswith('PARIS').values In [7]: # We create a vector of integers with the district # number of the corresponding station, or 0 # if the station is not in Paris. districts = np.zeros(len(paris), dtype=np.int32) districts[paris] = metro[4][paris].\ str.slice(6, 8).astype(np.int32) districts[~paris] = 0 ndistricts = districts.max() + 1 -
我们还提取所有地铁站的坐标：
In [8]: lon = metro[1] lat = metro[2] -
现在，让我们使用 OpenStreetMap 获取巴黎的地图。我们通过所有地铁站的极端纬度和经度坐标来指定地图的边界。我们使用轻量级的 Smopy 模块来生成地图：
In [9]: box = (lat[paris].min(), lon[paris].min(), lat[paris].max(), lon[paris].max()) m = smopy.Map(box, z=12) -
我们现在使用 SciPy 计算车站的 Voronoi 图。一个
Voronoi对象通过点的坐标创建。它包含了我们将用于显示的几个属性：In [10]: vor = spatial.Voronoi(np.c_[lat, lon]) -
我们创建了一个通用函数来显示 Voronoi 图。SciPy 已经实现了这样一个函数，但该函数没有考虑无限点。我们将使用的实现来自 Stack Overflow，地址是
stackoverflow.com/a/20678647/1595060。这个函数相对较长，我们不会在这里完全复制它。完整版本可以在书籍的 GitHub 仓库中找到。In [11]: def voronoi_finite_polygons_2d(vor, radius=None): """Reconstruct infinite Voronoi regions in a 2D diagram to finite regions.""" ... -
voronoi_finite_polygons_2d()函数返回一个区域列表和一个顶点列表。每个区域都是顶点索引的列表。所有顶点的坐标存储在vertices中。从这些结构中，我们可以创建一个单元格列表。每个单元格表示一个多边形，作为顶点坐标的数组。我们还使用smopy.Map实例的to_pixels()方法。此函数将纬度和经度的地理坐标转换为图像中的像素。In [12]: regions, vertices = \ voronoi_finite_polygons_2d(vor) cells = [m.to_pixels(vertices[region]) for region in regions] -
现在，我们计算每个多边形的颜色：
In [13]: cmap = plt.cm.Set3 # We generate colors for districts using # a color map. colors_districts = cmap( np.linspace(0., 1., ndistricts))[:,:3] # The color of every polygon, grey by default. colors = .25 * np.ones((len(districts), 3)) # We give each polygon in Paris the color of # its district. colors[paris] = colors_districts[districts[paris]] -
最后，我们使用
Map实例的show_mpl()方法来显示带有 Voronoi 图的地图：In [14]: ax = m.show_mpl() ax.add_collection( mpl.collections.PolyCollection(cells, facecolors=colors, edgecolors='k', alpha=0.35,))
它是如何工作的…
让我们给出欧几里得空间中 Voronoi 图的数学定义。如果*(x[i])是一个点集，那么这个点集的 Voronoi 图就是由以下定义的子集V[i]*（称为单元格或区域）的集合：
Voronoi 图的双重图是德劳内三角剖分。这个几何对象通过三角形覆盖了点集的凸包。
SciPy 使用Qhull（一个 C++计算几何库）来计算 Voronoi 图。
还有更多内容……
这里有更多参考资料：
-
维基百科上的 Voronoi 图，地址是
en.wikipedia.org/wiki/Voronoi_diagram -
维基百科上的德劳内三角剖分，地址是
en.wikipedia.org/wiki/Delaunay_triangulation -
scipy.spatial.voronoi文档可在docs.scipy.org/doc/scipy-dev/reference/generated/scipy.spatial.Voronoi.html找到 -
可在www.qhull.org获取 Qhull 库
另见
- 使用 Shapely 和 basemap 操作地理空间数据食谱
使用 Shapely 和 basemap 操作地理空间数据
在这个食谱中，我们将展示如何加载和显示 Shapefile 格式的地理数据。具体而言，我们将使用自然地球（www.naturalearthdata.com）的数据来显示非洲各国，并通过其人口和国内生产总值（GDP）进行着色。
Shapefile（en.wikipedia.org/wiki/Shapefile）是一种广泛使用的地理空间矢量数据格式，适用于 GIS 软件。它可以通过Fiona读取，Fiona 是GDAL/OGR（一个支持 GIS 文件格式的 C++库）的 Python 封装库。我们还将使用Shapely，这是一个用于处理二维几何图形的 Python 包，以及descartes，它用于在 matplotlib 中渲染 Shapely 图形。最后，我们将使用basemap绘制地图。
准备工作
你需要以下包：
-
GDAL/OGR 可以在www.gdal.org/ogr/找到
-
Fiona 可以在
toblerity.org/fiona/README.html找到 -
Shapely 可以在
toblerity.org/shapely/project.html找到 -
descartes 可以在
pypi.python.org/pypi/descartes找到 -
Basemap 可以在
matplotlib.org/basemap/找到
使用 Anaconda，你可以执行以下操作：
conda install gdal
conda install fiona
conda install basemap
Shapely 和 descartes 可以通过以下命令安装：
pip install shapely
pip install descartes
在 Windows 系统中，除了 descartes，你可以在 Chris Gohlke 的网页www.lfd.uci.edu/~gohlke/pyt…上找到所有这些包的二进制安装程序。
在其他系统中，你可以在项目网站上找到安装说明。GDAL/OGR 是 Fiona 所必需的 C++库。其他包是常规的 Python 包。
最后，你需要在本书的 GitHub 仓库中下载非洲数据集，链接为github.com/ipython-books/cookbook-data。数据来自 Natural Earth 网站，地址为www.naturalearthdata.com/downloads/1…。
如何操作……
-
让我们导入相关包：
In [1]: import numpy as np import matplotlib.pyplot as plt import matplotlib.collections as col from mpl_toolkits.basemap import Basemap import fiona import shapely.geometry as geom from descartes import PolygonPatch %matplotlib inline -
我们使用 Fiona 加载Shapefile数据集。这个数据集包含了全世界所有国家的边界。
In [2]: # Natural Earth data countries = fiona.open( "data/ne_10m_admin_0_countries.shp") -
我们选择非洲的国家：
In [3]: africa = [c for c in countries if c['properties']['CONTINENT'] == 'Africa'] -
现在，我们创建一个显示非洲大陆的底图：
In [4]: m = Basemap(llcrnrlon=-23.03, llcrnrlat=-37.72, urcrnrlon=55.20, urcrnrlat=40.58) -
让我们编写一个函数，将国家边界的地理坐标转换为地图坐标。这将使我们能够在底图中显示边界：
In [5]: def _convert(poly, m): if isinstance(poly, list): return [_convert(_, m) for _ in poly] elif isinstance(poly, tuple): return m(*poly) In [6]: for _ in africa: _['geometry']['coordinates'] = _convert( _['geometry']['coordinates'], m) -
下一步是创建来自 Fiona 加载的Shapefile数据集的 matplotlib
PatchCollection对象。我们将使用 Shapely 和 descartes 来完成：In [7]: def get_patch(shape, **kwargs): """Return a matplotlib PatchCollection from a geometry object loaded with fiona.""" # Simple polygon. if isinstance(shape, geom.Polygon): return col.PatchCollection( [PolygonPatch(shape, **kwargs)], match_original=True) # Collection of polygons. elif isinstance(shape, geom.MultiPolygon): return col.PatchCollection( [PolygonPatch(c, **kwargs) for c in shape], match_original=True) In [8]: def get_patches(shapes, fc=None, ec=None, **kwargs): """Return a list of matplotlib PatchCollection objects from a Shapefile dataset.""" # fc and ec are the face and edge colors of the # countries. We ensure these are lists of # colors, with one element per country. if not isinstance(fc, list): fc = [fc] * len(shapes) if not isinstance(ec, list): ec = [ec] * len(shapes) # We convert each polygon to a matplotlib # PatchCollection object. return [get_patch(geom.shape(s['geometry']), fc=fc_, ec=ec_, **kwargs) for s, fc_, ec_ in zip(shapes, fc, ec)] -
我们还定义了一个函数，根据Shapefile数据集中的特定字段获取国家的颜色。实际上，我们的数据集不仅包含国家边界，还包含每个国家的行政、经济和地理属性。在这里，我们将根据国家的人口和 GDP 选择颜色：
In [9]: def get_colors(field, cmap): """Return one color per country, depending on a specific field in the dataset.""" values = [country['properties'][field] for country in africa] values_max = max(values) return [cmap(v / values_max) for v in values] -
最后，我们显示地图。我们使用 basemap 显示海岸线，并使用我们的Shapefile数据集显示国家：
In [10]: # Display the countries color-coded with # their population. ax = plt.subplot(121) m.drawcoastlines() patches = get_patches(africa, fc=get_colors('POP_EST', plt.cm.Reds), ec='k') for p in patches: ax.add_collection(p) plt.title("Population") # Display the countries color-coded with # their population. ax = plt.subplot(122) m.drawcoastlines() patches = get_patches(africa, fc=get_colors('GDP_MD_EST', plt.cm.Blues), ec='k') for p in patches: ax.add_collection(p) plt.title("GDP")
另见
- 创建道路网络路线规划器的配方
创建一个道路网络的路线规划器
在这个实例中，我们基于前面几个实例中描述的技术，创建了一个简单的类似 GPS 的路线规划器。我们将从美国人口普查局获取加利福尼亚州的道路网络数据，用来在道路网络图中寻找最短路径。这使我们能够显示加利福尼亚州任意两地之间的道路路线。
准备工作
本实例需要 NetworkX 和 Smopy。在让 NetworkX 读取 Shapefile 数据集时，还需要 GDAL/OGR。更多信息可以参考前面的实例。
你还需要从书籍的 GitHub 仓库下载 Road 数据集：github.com/ipython-books/cookbook-data，并将其解压到当前目录。
注意
在写这篇文章时，NetworkX 对 Shapefile 的支持似乎与 Python 3.x 不兼容。因此，这个方法仅在 Python 2.x 下成功测试过。
如何实现…
-
让我们导入相关包：
In [1]: import networkx as nx import numpy as np import pandas as pd import json import smopy import matplotlib.pyplot as plt %matplotlib inline -
我们使用 NetworkX 加载数据（Shapefile 数据集）。这个数据集包含了加利福尼亚州主要道路的详细信息。NetworkX 的
read_shp()函数返回一个图，其中每个节点是一个地理位置，每条边包含连接两个节点的道路信息。数据来自美国人口普查局网站 www.census.gov/geo/maps-da…。In [2]: g = nx.read_shp("data/tl_2013_06_prisecroads.shp") -
这个图不一定是连通的，但为了计算最短路径，我们需要一个连通图。在这里，我们使用
connected_component_subgraphs()函数获取最大的连通子图：In [3]: sgs = list(nx.connected_component_subgraphs( g.to_undirected())) largest = np.argmax([len(sg) for sg in sgs]) sg = sgs[largest] len(sg) Out[3]: 464 -
我们定义了两个位置（包含纬度和经度），并找出这两个位置之间的最短路径：
In [4]: pos0 = (36.6026, -121.9026) pos1 = (34.0569, -118.2427) -
图中的每条边包含道路的详细信息，包括沿道路的点的列表。我们首先创建一个函数，返回图中任意一条边的坐标数组：
In [5]: def get_path(n0, n1): """If n0 and n1 are connected nodes in the graph, this function returns an array of point coordinates along the road linking these two nodes.""" return np.array(json.loads( sg[n0][n1]['Json'])['coordinates']) -
我们可以通过道路路径来计算其长度。我们首先需要定义一个函数，计算地理坐标系中任意两点之间的距离。这个函数可以在 Stack Overflow 上找到 (
stackoverflow.com/questions/8858838/need-help-calculating-geographical-distance)：In [6]: EARTH_R = 6372.8 def geocalc(lat0, lon0, lat1, lon1): """Return the distance (in km) between two points in geographical coordinates.""" lat0 = np.radians(lat0) lon0 = np.radians(lon0) lat1 = np.radians(lat1) lon1 = np.radians(lon1) dlon = lon0 - lon1 y = np.sqrt( (np.cos(lat1)*np.sin(dlon))**2 +(np.cos(lat0)*np.sin(lat1) -np.sin(lat0)*np.cos(lat1)* \ np.cos(dlon))**2) x = np.sin(lat0)*np.sin(lat1) + \ np.cos(lat0)*np.cos(lat1)*np.cos(dlon) c = np.arctan2(y, x) return EARTH_R * c -
现在，我们定义一个计算路径长度的函数：
In [7]: def get_path_length(path): return np.sum(geocalc( path[1:,0], path[1:,1], path[:-1,0], path[:-1,1])) -
现在，我们通过计算任何两个连接节点之间的距离来更新我们的图。我们将这个信息添加到边的
distance属性中：In [8]: # Compute the length of the road segments. for n0, n1 in sg.edges_iter(): path = get_path(n0, n1) distance = get_path_length(path) sg.edge[n0][n1]['distance'] = distance -
在我们能够找到图中的最短路径之前，最后一步是找出图中最接近两个请求位置的两个节点：
In [9]: nodes = np.array(sg.nodes()) # Get the closest nodes in the graph. pos0_i = np.argmin(np.sum( (nodes[:,::-1] - pos0)**2, axis=1)) pos1_i = np.argmin(np.sum( (nodes[:,::-1] - pos1)**2, axis=1)) -
现在，我们使用 NetworkX 的
shortest_path()函数来计算两个位置之间的最短路径。我们指定每条边的权重为它们之间道路的长度：In [10]: # Compute the shortest path. path = nx.shortest_path(sg, source=tuple(nodes[pos0_i]), target=tuple(nodes[pos1_i]), weight='distance') len(path) Out[10]: 19 -
行程已计算完成。
path变量包含了构成我们两点之间最短路径的边的列表。现在，我们可以通过 pandas 获取有关行程的信息。数据集中有几个感兴趣的字段，包括道路的名称和类型（州道、州际公路等）：In [11]: roads = pd.DataFrame([ sg.edge[path[i]][path[i + 1]] for i in range(len(path)-1)], columns=['FULLNAME', 'MTFCC', 'RTTYP', 'distance']) roads Out[11]: FULLNAME MTFCC RTTYP distance 0 State Rte 1 S1200 S 100.657768 1 State Rte 1 S1200 S 33.419581 ... 16 Hollywood Fwy S1200 M 14.087627 17 Hollywood Fwy S1200 M 0.010107这是此行程的总长度：
In [12]: roads['distance'].sum() Out[12]: 508.66421585288725 -
最后，让我们在地图上展示该行程。我们首先通过 Smopy 获取地图：
In [13]: map = smopy.Map(pos0, pos1, z=7, margin=.1) -
我们的路径包含图中的连接节点。每两个节点之间的边由一系列点（构成道路的一部分）来描述。因此，我们需要定义一个函数，将路径中每个边上的位置串联起来。我们必须按照正确的顺序连接这些位置。我们选择的顺序是基于这样一个事实：一条边的最后一个点需要接近下一条边的第一个点：
In [14]: def get_full_path(path): """Return the positions along a path.""" p_list = [] curp = None for i in range(len(path)-1): p = get_path(path[i], path[i+1]) if curp is None: curp = p if np.sum((p[0]-curp)**2) > \ np.sum((p[-1]-curp)**2): p = p[::-1,:] p_list.append(p) curp = p[-1] return np.vstack(p_list) -
我们将路径转换为像素，以便在 Smopy 地图上显示：
In [15]: linepath = get_full_path(path) x, y = map.to_pixels(linepath[:,1], linepath[:,0]) -
最后，让我们展示地图，标出我们的两个位置以及它们之间计算出的行程：
In [16]: map.show_mpl() # Plot the itinerary. plt.plot(x, y, '-k', lw=1.5) # Mark our two positions. plt.plot(x[0], y[0], 'ob', ms=10) plt.plot(x[-1], y[-1], 'or', ms=10)
它是如何工作的…
我们使用 NetworkX 的 shortest_path() 函数计算了最短路径。这里，这个函数使用了Dijkstra 算法。这个算法有着广泛的应用，例如网络路由协议。
计算两个点之间的地理距离有多种方法。在这里，我们使用了一个相对精确的公式：正弦距离（也叫做大圆距离），该公式假设地球是一个完美的球体。我们也可以使用一个更简单的公式，因为在一条路上两个连续点之间的距离很小。
还有更多内容…
您可以在以下参考资料中找到关于最短路径问题和 Dijkstra 算法的更多信息：
-
维基百科上的最短路径，
en.wikipedia.org/wiki/Shortest_path_problem -
Dijkstra 算法，描述见
en.wikipedia.org/wiki/Dijkstra%27s_algorithm
这里有一些关于地理距离的参考资料：
-
维基百科上的地理距离，
en.wikipedia.org/wiki/Geographical_distance -
维基百科上的大圆，
en.wikipedia.org/wiki/Great_circle -
维基百科上的大圆距离，
en.wikipedia.org/wiki/Great-circle_distance
第十五章：符号和数值数学
在本章中，我们将涵盖以下主题：
-
使用 SymPy 进行符号计算
-
解方程和不等式
-
分析实值函数
-
计算精确概率和操作随机变量
-
用 SymPy 进行一点数论
-
从真值表中找到布尔命题公式
-
分析非线性微分系统 - Lotka-Volterra（捕食者-猎物）方程
-
开始使用 Sage
介绍
在本章中，我们将介绍SymPy，这是一个用于符号数学的 Python 库。虽然本书大部分内容涉及数值方法，但在这里我们将看到符号计算更适合的示例。
SymPy 对符号计算就像 NumPy 对数值计算一样重要。例如，SymPy 可以帮助我们在运行模拟之前分析数学模型。
尽管 SymPy 相当强大，但与其他计算机代数系统相比有点慢。主要原因是 SymPy 是用纯 Python 编写的。一个更快更强大的数学系统是Sage（也请参阅本章中的Getting started with Sage示例）。Sage 是一个庞大的独立程序，有许多大型依赖项（包括 SymPy！），并且在撰写本文时仅使用 Python 2。它主要用于交互式使用。Sage 包括类似 IPython 的笔记本。
LaTeX
LaTeX是一种广泛用于编写出版质量数学方程式的文档标记语言。使用 LaTeX 编写的方程式可以在浏览器中使用MathJax JavaScript 库显示。SymPy 使用此系统在 IPython 笔记本中显示方程式。
LaTeX 方程式也可以在 matplotlib 中使用。在这种情况下，建议在本地计算机上安装 LaTeX。
这里有一些参考资料：
-
Wikipedia 上的 LaTeX，网址为
en.wikipedia.org/wiki/LaTeX -
MathJax，可在www.mathjax.org找到
-
matplotlib 中的 LaTeX，描述在
matplotlib.org/users/usetex.html -
用于显示 SymPy 方程式的文档，网址为
docs.sympy.org/latest/tutorial/printing.html -
要在计算机上安装 LaTeX，请参考
latex-project.org/ftp.html
使用 SymPy 进行符号计算
在这个示例中，我们将简要介绍使用 SymPy 进行符号计算。我们将在接下来的示例中看到 SymPy 的更高级功能。
准备工作
SymPy 是一个纯 Python 包，没有其他依赖项，因此非常容易安装。使用 Anaconda，您可以在终端中键入conda install sympy。在 Windows 上，您可以使用 Chris Gohlke 的软件包(www.lfd.uci.edu/~gohlke/pyt…)。最后，您可以使用pip install sympy命令。
如何做...
SymPy 可以在 Python 模块中使用，或者在 IPython 中交互式使用。在笔记本中，所有数学表达式都通过 LaTeX 显示，得益于 MathJax JavaScript 库。
这里是 SymPy 的介绍：
-
首先，我们导入 SymPy 并在 IPython 笔记本中启用 LaTeX 打印：
In [1]: from sympy import * init_printing() -
要处理符号变量，我们首先需要声明它们：
In [2]: var('x y') Out[2]: (x, y) -
var()函数用于创建符号并将其注入到命名空间中。这个函数仅应在交互模式下使用。在 Python 模块中，最好使用symbols()函数来返回符号：In [3]: x, y = symbols('x y') -
我们可以用这些符号创建数学表达式：
In [4]: expr1 = (x + 1)**2 expr2 = x**2 + 2*x + 1 -
这些表达式相等吗？
In [5]: expr1 == expr2 Out[5]: False -
这些表达式在数学上相等，但在语法上并不相同。为了测试它们是否在数学上相等，我们可以让 SymPy 代数化地简化差异：
In [6]: simplify(expr1-expr2) Out[6]: 0 -
处理符号表达式时，一个非常常见的操作是通过符号表达式的
subs()方法将一个符号替换为另一个符号、表达式或数字：在 SymPy 表达式中的替代
-
一个有理数不能简单地写作
1/2，因为这个 Python 表达式会计算出 0。一个解决方法是将数字1转换为 SymPy 整数对象，例如使用S()函数：In [9]: expr1.subs(x, S(1)/2) Out[9]: 9/4 -
精确表示的数字可以通过
evalf进行数值求值：In [10]: _.evalf() Out[10]: 2.25000000000000 -
我们可以通过
lambdify()函数轻松地从 SymPy 符号表达式创建 Python 函数。生成的函数可以特别地在 NumPy 数组上进行求值。这在我们需要从符号世界转到数值世界时非常方便：In [11]: f = lambdify(x, expr1) In [12]: import numpy as np f(np.linspace(-2., 2., 5)) Out[12]: array([ 1., 0., 1., 4., 9.])
它是如何工作的...
SymPy 的一个核心思想是使用标准的 Python 语法来操作精确的表达式。虽然这非常方便和自然，但也有一些注意事项。像x这样的符号，代表数学变量，在被实例化之前不能在 Python 中使用（否则，解释器会抛出NameError异常）。这与大多数其他计算机代数系统不同。因此，SymPy 提供了提前声明符号变量的方法。
另一个例子是整数除法；由于1/2会计算为0（在 Python 2 中），SymPy 无法知道用户原本打算写一个分数。我们需要先将数值整数1转换为符号整数1，然后再除以2。
此外，Python 的等式指的是语法树之间的相等，而不是数学表达式之间的相等。
另见
-
求解方程和不等式食谱
-
Sage 入门食谱
求解方程和不等式
SymPy 提供了多种方法来求解线性和非线性方程及方程组。当然，这些函数并不总能成功找到封闭形式的精确解。在这种情况下，我们可以退回到数值求解器并得到近似解。
准备就绪
我们首先需要导入 SymPy。我们还需要在笔记本中初始化美观打印（参见本章的第一个示例）。
如何实现...
-
定义一些符号：
In [2]: var('x y z a') Out[2]: (x, y, z, a) -
我们使用
solve()函数来解方程（右侧默认是0）：In [3]: solve(x**2 - a, x) Out[3]: [-sqrt(a), sqrt(a)] -
我们还可以解不等式。在这里，我们需要使用
solve_univariate_inequality()函数来解这个实数域中的单变量不等式：In [4]: x = Symbol('x') solve_univariate_inequality(x**2 > 4, x) Out[4]: Or(x < -2, x > 2) -
solve()函数也接受方程组（这里是一个线性系统）：In [5]: solve([x + 2*y + 1, x - 3*y - 2], x, y) Out[5]: {x: 1/5, y: -3/5} -
非线性系统也可以处理：
In [6]: solve([x**2 + y**2 - 1, x**2 - y**2 - S(1)/2], x, y) Out[6]: [(-sqrt(3)/2, -1/2), (-sqrt(3)/2, 1/2), (sqrt(3)/2, -1/2), (sqrt(3)/2, 1/2)] -
奇异线性系统也可以求解（这里有无限多个解，因为两个方程是共线的）：
In [7]: solve([x + 2*y + 1, -x - 2*y - 1], x, y) Out[7]: {x: -2*y - 1} -
现在，让我们使用包含符号变量的矩阵来求解线性系统：
In [8]: var('a b c d u v') Out[8]: (a, b, c, d, u, v) -
我们创建了增广矩阵，这是系统矩阵与线性系数和右侧向量的水平拼接。这个矩阵对应于以下系统 x,y: ax+by=u, cx+dy=v：
In [9]: M = Matrix([[a, b, u], [c, d, v]]); M Out[9]: Matrix([[a, b, u], [c, d, v]]) In [10]: solve_linear_system(M, x, y) Out[10]: {x: (-b*v + d*u)/(a*d - b*c), y: ( a*v - c*u)/(a*d - b*c)} -
这个系统需要是非奇异的，才能有唯一解，这等同于说系统矩阵的行列式需要非零（否则前述分数的分母将为零）：
In [11]: det(M[:2,:2]) Out[11]: a*d - b*c
还有更多内容...
SymPy 中的矩阵支持非常丰富；我们可以执行大量的操作和分解（请参阅docs.sympy.org/latest/modules/matrices/matrices.html的参考指南）。
这里有更多关于线性代数的参考资料：
-
维基百科上的线性代数，
en.wikipedia.org/wiki/Linear_algebra#Further_reading -
维基教科书上的线性代数，
en.wikibooks.org/wiki/Linear_Algebra
分析实值函数
SymPy 包含了丰富的微积分工具箱，用于分析实值函数：极限、幂级数、导数、积分、傅里叶变换等等。在这个示例中，我们将展示这些功能的基本用法。
做好准备
我们首先需要导入 SymPy。我们还需要在笔记本中初始化美观打印（参见本章的第一个示例）。
如何实现...
-
定义一些符号和一个函数（它只是一个依赖于
x的表达式）：In [1]: var('x z') Out[1]: (x, z) In [2]: f = 1/(1+x**2) -
让我们在
1处评估这个函数：In [3]: f.subs(x, 1) Out[3]: 1/2 -
我们可以计算这个函数的导数：
In [4]: diff(f, x) Out[4]: -2*x/(x**2 + 1)**2 -
f的极限是无限大吗？（注意符号的双oo（无穷大）表示）：In [5]: limit(f, x, oo) Out[5]: 0 -
下面是如何计算泰勒级数（这里是围绕
0，阶数为9）。大 O 可以通过removeO()方法去除。In [6]: series(f, x0=0, n=9) Out[6]: 1 - x**2 + x**4 - x**6 + x**8 + O(x**9) -
我们可以计算定积分（这里是对整个实数轴的积分）：
In [7]: integrate(f, (x, -oo, oo)) Out[7]: pi -
SymPy 还可以计算不定积分：
In [8]: integrate(f, x) Out[8]: atan(x) -
最后，让我们计算
f的傅里叶变换：In [9]: fourier_transform(f, x, z) Out[9]: pi*exp(-2*pi*z)
还有更多内容...
SymPy 还包括许多其他的积分变换，除了傅里叶变换（docs.sympy.org/dev/modules/integrals/integrals.html）。然而，SymPy 并不总是能够找到闭式解。
这里有一些关于实分析和微积分的参考书目：
-
维基百科上的实分析内容，访问
en.wikipedia.org/wiki/Real_analysis#Bibliography -
维基书上的微积分内容，访问
en.wikibooks.org/wiki/Calculus
计算精确概率并操作随机变量
SymPy 包括一个名为stats的模块，允许我们创建和操作随机变量。当我们处理概率或统计模型时，这非常有用；我们可以计算符号期望值、方差、概率和随机变量的密度。
如何操作...
-
我们来导入 SymPy 和 stats 模块：
In [1]: from sympy import * from sympy.stats import * init_printing() -
让我们掷两个骰子，
X和Y，每个都有六个面：In [2]: X, Y = Die('X', 6), Die('Y', 6) -
我们可以计算由等式（使用
Eq运算符）或不等式定义的概率：In [3]: P(Eq(X, 3)) Out[3]: 1/6 In [4]: P(X>3) Out[4]: 1/2 -
条件也可以涉及多个随机变量：
In [5]: P(X>Y) Out[5]: 5/12 -
我们可以计算条件概率：
In [6]: P(X+Y>6, X<5) Out[6]: 5/12 -
我们也可以处理任意的离散或连续随机变量：
In [7]: Z = Normal('Z', 0, 1) # Gaussian variable In [8]: P(Z>pi) Out[8]: -erf(sqrt(2)*pi/2)/2 + 1/2 -
我们可以计算期望值和方差：
In [9]: E(Z**2), variance(Z**2) Out[9]: (1, 2) -
我们也可以计算密度：
In [10]: f = density(Z) In [11]: var('x') f(x) Out[11]: sqrt(2)*exp(-x**2/2)/(2*sqrt(pi)) -
我们可以绘制这些密度：
In [12]: %matplotlib inline plot(f(x), (x, -6, 6))高斯密度
它是如何工作的...
SymPy 的stats模块包含许多函数，用于定义具有经典分布（如二项分布、指数分布等）的随机变量，无论是离散的还是连续的。它通过利用 SymPy 强大的积分算法来计算概率分布的积分，从而精确计算概率量。例如，就是：
请注意，等式条件是使用Eq运算符而非更常见的 Python 语法==来表示的。这是 SymPy 的一个通用特性；==表示 Python 变量之间的相等，而Eq表示符号表达式之间的数学运算。
一点关于数论的内容与 SymPy
SymPy 包含许多与数论相关的例程：获取质数、整数分解等等。我们将在这里展示几个例子。
准备就绪
要在 matplotlib 中使用 LaTeX 显示图例，你需要在计算机上安装 LaTeX（请参见本章的简介）。
如何操作...
-
我们来导入 SymPy 和数论包：
In [1]: from sympy import * init_printing() In [2]: import sympy.ntheory as nt -
我们可以测试一个数字是否是质数：
In [3]: nt.isprime(2011) Out[3]: True -
我们可以找出给定数字后的下一个质数：
In [4]: nt.nextprime(2011) Out[4]: 2017 -
第 1000 个质数是什么？
In [5]: nt.prime(1000) Out[5]: 7919 -
小于 2011 的质数有多少个？
In [6]: nt.primepi(2011) Out[6]: 305 -
我们可以绘制
，素数计数函数（小于或等于某个数字 x 的素数的数量）。著名的素数定理指出，这个函数在渐近意义上等价于 x/log(x)。这个表达式大致量化了素数在所有整数中的分布情况：
In [7]: import numpy as np import matplotlib.pyplot as plt %matplotlib inline x = np.arange(2, 10000) plt.plot(x, map(nt.primepi, x), '-k', label='$\pi(x)$') plt.plot(x, x / np.log(x), '--k', label='$x/\log(x)$') plt.legend(loc=2)素数分布
-
让我们计算一个数字的整数因式分解：
In [8]: nt.factorint(1998) Out[8]: {2: 1, 3: 3, 37: 1} In [9]: 2 * 3**3 * 37 Out[9]: 1998 -
最后，一个小问题。一个懒惰的数学家在数他的珠子。当它们排列成三行时，最后一列有一颗珠子。当它们排列成四行时，最后一列有两颗珠子，排列成五行时有三颗珠子。那么珠子总共有多少颗？（提示：懒惰的数学家少于 100 颗珠子。）
使用中国剩余定理计数珠子
中国剩余定理给我们提供了答案：
In [10]: from sympy.ntheory.modular import solve_congruence In [11]: solve_congruence((1, 3), (2, 4), (3, 5)) Out[11]: (58, 60)有无限多的解：58 加上任何 60 的倍数。由于珠子总数少于 100 颗，58 就是正确答案。
它是如何工作的…
SymPy 包含许多与数论相关的函数。在这里，我们使用中国剩余定理来求解以下算术方程组的解：
中国剩余定理
三重条是模同余的符号。这里，它表示 m[i] 整除 a[i]-n。换句话说，n 和 a[i] 相等，直到 m[i] 的倍数。处理同余时，当涉及周期性尺度时非常方便。例如，12 小时制的时钟操作是按模 12 运算的。数字 11 和 23 在模 12 下是等价的（它们在时钟上表示相同的时刻），因为它们的差是 12 的倍数。
在这个例子中，必须满足三个同余：珠子数在除以 3 时的余数是 1（表示这种排列中多了一颗珠子），除以 4 时余数是 2，除以 5 时余数是 3。使用 SymPy，我们只需在solve_congruence()函数中指定这些值，就能得到解。
该定理指出，当 m[i] 彼此互质时（它们之间的任意两个数都是互质的），解是存在的。所有解在 m[i] 的乘积下是同余的。这个数论中的基本定理有许多应用，特别是在密码学中。
还有更多…
这里有一些关于数论的教材：
-
本科水平：《初等数论》，Gareth A. Jones，Josephine M. Jones，Springer，(1998)
-
研究生水平：《现代数论的经典介绍》，Kenneth Ireland，Michael Rosen，Springer，(1982)
这里有一些参考资料：
-
SymPy 的数论模块文档，访问地址：
docs.sympy.org/dev/modules/ntheory.html -
中国剩余定理的应用，见于
mathoverflow.net/questions/10014/applications-of-the-chinese-remainder-theorem
从真值表中找到一个布尔命题公式
SymPy 中的逻辑模块让我们操作复杂的布尔表达式，也就是命题公式。
这个例子将展示这个模块的一个实际应用。假设在一个程序中，我们需要根据三个布尔变量写一个复杂的if语句。我们可以考虑八种可能的情况（真、真与假、依此类推），并评估每种情况下应该得到什么结果。SymPy 提供了一个函数，可以生成一个满足我们真值表的紧凑逻辑表达式。
如何操作...
-
让我们导入 SymPy：
In [1]: from sympy import * init_printing() -
让我们定义一些符号：
In [2]: var('x y z') -
我们可以使用符号和几个运算符来定义命题公式：
In [3]: P = x & (y | ~z); P Out[3]: And(Or(Not(z), y), x) -
我们可以使用
subs()来对实际的布尔值进行公式求值：In [4]: P.subs({x: True, y: False, z: True}) Out[4]: False -
现在，我们希望根据
x、y和z找到一个命题公式，以下是其真值表：一张真值表
-
让我们列出所有我们希望求值为
True的组合，以及那些结果无关紧要的组合：In [6]: minterms = [[1,0,1], [1,0,0], [0,0,0]] dontcare = [[1,1,1], [1,1,0]] -
现在，我们使用
SOPform()函数来推导出一个合适的公式：In [7]: Q = SOPform(['x', 'y', 'z'], minterms, dontcare); Q Out[7]: Or(And(Not(y), Not(z)), x) -
让我们测试一下这个命题是否有效：
In [8]: Q.subs({x: True, y: False, z: False}), Q.subs({x: False, y: True, z: True}) Out[8]: (True, False)
它是如何工作的...
SOPform()函数生成一个与真值表对应的完整表达式，并使用Quine-McCluskey 算法简化它。它返回最小的积和形式（或合取的析取）。类似地，POSform()函数返回一个和的积。
给定的真值表可能出现在这种情况下：假设我们希望在文件不存在时写入文件（z），或者用户希望强制写入（x）。此外，用户可以阻止写入（y）。当文件需要被写入时，表达式的值为True。得到的 SOP 公式在我们首先显式禁止x和y时有效（强制和禁止写入同时发生是禁止的）。
还有更多内容...
这里有一些参考资料：
分析非线性微分系统 – Lotka-Volterra（捕食-被捕食）方程
在这里，我们将对著名的非线性微分方程系统——洛特卡-沃尔特拉方程进行简要的分析，这也被称为捕食者-猎物方程。这些方程是描述两种相互作用的种群（例如鲨鱼和沙丁鱼）演化的一级微分方程，其中捕食者捕食猎物。这个例子展示了如何使用 SymPy 获得关于固定点及其稳定性的精确表达式和结果。
准备好
对于这个公式，建议了解线性和非线性微分方程的基础知识。
如何做到...
-
让我们创建一些符号：
In [1]: from sympy import * init_printing() In [2]: var('x y') var('a b c d', positive=True) Out[2]: (a, b, c, d) -
变量
x和y分别代表猎物和捕食者的种群。参数a、b、c和d是严格为正的参数（在这个公式的*如何实现...*部分中更精确地描述）。方程为：洛特卡-沃尔特拉方程
In [3]: f = x * (a - b*y) g = -y * (c - d*x) -
让我们找到系统的固定点（解f(x,y) = g(x,y) = 0）。我们将它们称为*(x[0], y[0])和(x[1], y[1])*：
In [4]: solve([f, g], (x, y)) Out[4]: [(0, 0), (c/d, a/b)] In [5]: (x0, y0), (x1, y1) = _ -
让我们用两个方程来写出二维向量：
In [6]: M = Matrix((f, g)); M Out[6]: Matrix([[ x*(a - b*y)], [-y*(c - d*x)]]) -
现在，我们可以计算系统的雅可比矩阵，作为
(x, y)的函数：In [7]: J = M.jacobian((x, y)); J Out[7]: Matrix([ [a - b*y, -b*x], [ d*y, -c + d*x]]) -
让我们通过观察雅可比矩阵在这个点的特征值来研究第一个固定点的稳定性。第一个固定点对应的是种群灭绝：
In [8]: M0 = J.subs(x, x0).subs(y, y0); M0 Out[8]: Matrix([a, 0], [0, -c]]) In [9]: M0.eigenvals() Out[9]: {a: 1, -c: 1}参数
a和c是严格为正的，因此特征值是实数且符号相反，这个固定点是鞍点。由于这个点不稳定，因此在这个模型中两种种群都灭绝的可能性很小。 -
现在让我们考虑第二个固定点：
In [10]: M1 = J.subs(x, x1).subs(y, y1); M1 Out[10]: Matrix([[ 0, -b*c/d], [a*d/b, 0]]) In [11]: M1.eigenvals() Out[11]: {-I*sqrt(a)*sqrt(c): 1, I*sqrt(a)*sqrt(c): 1}特征值是纯虚数；因此，这个固定点不是双曲的。因此，我们不能从这个线性分析中得出关于该固定点周围系统定性行为的结论。然而，我们可以通过其他方法证明，在这个点附近会发生振荡。
如何实现...
洛特卡-沃尔特拉方程模拟了捕食者和猎物种群的增长，考虑了它们的相互作用。在第一个方程中，ax项表示猎物的指数增长，-bxy表示捕食者导致的死亡。同样，在第二个方程中，-yc表示捕食者的自然死亡，dxy表示捕食者因捕食更多猎物而增长。
要找到系统的平衡点，我们需要找到* x, y 的值，使得 dx/dt = dy/dt = 0*，即f(x, y) = g(x, y) = 0，这样变量就不再变化了。在这里，我们能够通过solve()函数获得这些平衡点的解析值。
为了分析它们的稳定性，我们需要对非线性方程进行线性分析，通过在这些平衡点计算雅可比矩阵。这个矩阵代表了线性化系统，它的特征值告诉我们关于平衡点附近系统的稳定性。哈特曼–格罗曼定理指出，如果平衡点是双曲的（意味着矩阵的特征值没有实部为零的情况），那么原始系统的行为在定性上与线性化系统在平衡点附近的行为是匹配的。在这里，第一个平衡点是双曲的，因为a, c > 0，而第二个则不是。
在这里，我们能够计算出平衡点处雅可比矩阵及其特征值的符号表达式。
还有更多...
即使一个差分系统无法解析求解（如这里所示），数学分析仍然可以为我们提供关于系统解的行为的定性信息。当我们关心定性结果时，纯粹的数值分析并不总是相关的，因为数值误差和近似可能导致关于系统行为的错误结论。
这里有一些参考资料：
-
SymPy 中的矩阵文档，网址：
docs.sympy.org/dev/modules/matrices/matrices.html -
维基百科上的动态系统，
en.wikipedia.org/wiki/Dynamical_system -
Scholarpedia 上的平衡点，www.scholarpedia.org/article/Equ…
-
维基百科上的分岔理论，
en.wikipedia.org/wiki/Bifurcation_theory -
维基百科上的混沌理论，
en.wikipedia.org/wiki/Chaos_theory -
关于动态系统的进一步阅读，网址：
en.wikipedia.org/wiki/Dynamical_system#Further_reading
开始使用 Sage
Sage (www.sagemath.org) 是一个基于 Python 的独立数学软件。它是商业产品如 Mathematica、Maple 或 MATLAB 的开源替代品。Sage 提供了一个统一的接口，连接了许多开源数学库。这些库包括 SciPy、SymPy、NetworkX 和其他 Python 科学包，也包括非 Python 库，如 ATLAS、BLAS、GSL、LAPACK、Singular 等。
在这个教程中，我们将简要介绍 Sage。
准备工作
你可以选择：
-
在本地计算机上安装 Sage (www.sagemath.org/doc/install…)
-
在云端远程创建 Sage 笔记本 (
cloud.sagemath.com/)
由于依赖于如此多的库，Sage 很庞大且难以从源代码编译。除了 Windows 系统外，其他大多数系统都有现成的二进制文件，而在 Windows 上通常需要使用 VirtualBox（一种虚拟化解决方案：www.virtualbox.org）。
另外，你也可以通过在云端运行 IPython 笔记本，在浏览器中使用 Sage。
请注意，Sage 在写作时与 Python 3 不兼容。
通常，Sage 是通过内置的笔记本进行交互式使用（它类似于 IPython 笔记本）。如果你想在 Python 程序中使用 Sage（也就是从 Python 导入 Sage），你需要运行 Sage 的内置 Python 解释器（www.sagemath.org/doc/faq/faq…）。
如何做...
在这里，我们将创建一个新的 Sage 笔记本，并介绍最基本的功能：
-
Sage 接受数学表达式，就像我们预期的那样：
sage: 3*4 12 -
由于基于 Python，Sage 的语法几乎与 Python 相同，但也有一些差异。例如，幂指数使用的是更经典的
^符号：sage: 2³ 8 -
就像在 SymPy 中一样，符号变量需要事先使用
var()函数声明。然而，x变量总是预定义的。在这里，我们定义一个新的数学函数：sage: f=1-sin(x)² -
让我们简化
f的表达式：sage: f.simplify_trig() cos(x)² -
让我们在给定的点上评估
f：sage: f(x=pi) 1 -
函数可以进行微分和积分：
sage: f.diff(x) -2*cos(x)*sin(x) sage: f.integrate(x) 1/2*x + 1/4*sin(2*x) -
除了符号计算，Sage 还支持数值计算：
sage: find_root(f-x, 0, 2) 0.6417143708729723 -
Sage 还具有丰富的绘图功能（包括交互式绘图控件）：
sage: f.plot((x, -2*pi, 2*pi))
还有更多内容...
这个（也）简短的教程无法充分展示 Sage 提供的广泛功能。Sage 涉及数学的许多方面：代数、组合数学、数值数学、数论、微积分、几何、图论等。以下是一些参考资料：
-
关于 Sage 的深入教程可以在 www.sagemath.org/doc/tutoria… 查阅
-
Sage 的参考手册可以在 www.sagemath.org/doc/referen… 查阅
-
关于 Sage 的视频，可以在 www.sagemath.org/help-video.… 查阅
另见
- 深入学习符号计算与 SymPy 课程
