| 姓名 | 学号 | 班级 | 选题 | 论述 | 结论 | 总分 |
| ---- | ---- | ---- |
| 余康 | 2014301020117 | 弘毅班 | 随机系统 | | | | |
不同类型随机行走过程的python模拟
<br >
余康 2014301020117 14级物理学弘毅班
<br >
前言
这篇文章是作者（我）在阅读过《计算物理》教材第七章“随机系统”之后的作品，由于时间与知识水平的限制，所以无法涵盖第七章的全部内容。基于以上的考虑，作者选择了相对较为基础的内容作为自己的课题，即“随机行走的过程的python模拟”。尽管如此，但是作者还是尽力地完善了自己的作品，相当于对书中第七章第一二节的内容做了自己的补充。
摘要
从布朗运动我们开始接触随机行走，布朗运动是由于运动是由于颗粒受液体分子碰撞的不平衡力引起的。随机行走可以理解为布朗运动的理想状态，也可以称为无规则行走，或者是形象的解释为醉汉的行走。不止于布朗运动，随机行走过程在我们生活中的实例与应用是极其广泛的，比如实例有所有的扩散现象，应用于互联网的链接分析，金融的股票分析和高分子的构象（自回避行走）等领域。随机行走过程的模拟可以通过随机数生成器来实现，从而我们可以通过数值方法验证解析方法得出随机系统的概率描述（例如平均数，方差等）并且直观地认识随机系统的规律。对于一维与二维随机行走的模拟较为简单，但是对于三维的随机行走（全方向，不限于格点）的模拟，就涉及如何在三维球面上随机取点（均匀分布），我采取了三维球面上的Marsaglia方法。
关键词： 随机行走 模拟 随机数 概率描述 Marsaglia方法
Ⅰ 介绍
随机行走是一个数学对象，它描述了由一系列随机步骤组成的路径。例如，跟踪液体或气体中分子的运动路径，搜索猎物的觅食动物的运动路径，超弦的行为，一支股票价格的波动和一个赌徒的财务情况都可以通过随机行走模型近似，尽管在现实中他们可能无法实现真正的随机。这些例子说明了随机行走在许多科学领域的应用， 包括生态学、心理学、计算机科学、物理学、化学和生物学，以及经济学。随机行走可以解释在这些领域的许多过程中所观察到的行为，从而作为随机活动的基本模型。谈到更数学的应用，π的值可以采用基于随机行走的编程方法数值逼近（Monte Carlo method，蒙特卡罗算法）。随机行走是由卡尔·皮尔逊在1905年第一次提出的。
人们对于随机行走有各种各样的不同方面的研究兴趣。这个词本身通常指的是一种特殊的马尔可夫链或马尔可夫过程。但是许多具有特定性质的随时间变化的过程也被称为随机行走。随机行走（无论是否马尔科夫）也能够发生在各种空间：大家通常研究的有图（图论），整数（数论），实线、平面、更高维的向量空间，曲面、更高维的黎曼流形，有限群、无限群、李群（群论）等。时间参数也可以被操控。在最简单的情形下，我们研究的是离散时间序列的随机行走。当然，我们也可以定义任意时间序列的随机行走，在这种情形下，所有时刻的位置都有定义。随机行走的实例包括莱维飞行，扩散模型（比如布朗运动）。
随机行走是马尔可夫过程讨论的一个基本主题.。它们的数学研究已十分广泛。随机行走的几个特性，包括分布（dispersal distributions），首次超越（first-passage），碰撞次数（ hitting times），会遇率（encounter rates），循环（recurrence），瞬态（transience），已被引入来量化它们的行为。
-
随机行走所满足的扩散定律：
其中D为扩散常数。
以下是四幅图是维基百科上对二维与三维随机行走的模拟结果，也是我在编程过程中力求达到的目标。
二维：
三维：
Ⅱ 正文
一维随机行走的模拟
一维的随机行走，是最简单的随机行走模型，但是它的研究方法以及某些结论却并不局限于一维的情况，可以轻易地平行推广到更高维的情形，所以一维随机行走的讨论是本文的基础。
定步长，等概率
考虑最简单的情形，即一维定步长随机行走，也就是左右等可能随机行走，我们把模型简化为：一个醉汉在一条直线上行走，每次只能走一个单位长度，并且他向每一边走的概率都是完全相等的（1/2）。我采用random.random()函数生成一个从0到1的随机数作为醉汉每一步行走方向的标尺，如果这个随机数大于0.5，则醉汉下一步会向右走一个单位长度；否则，醉汉下一步会向左走一个单位长度。
借助模拟，我们可以观察醉汉行走一百步之后的位移，但是由于每一次醉汉的行走都是随机的，有意义的只是模拟次数足够多之后所表现出来的统计规律，所以我们可以假想有一个醉汉系综，实际上我们通过循环结构即可实现，他们都行走了一百步之后，对他们与原点的位移取平均值，才能验证我们对这个随机系统（醉汉系综）的性质的预测。
接下来我将展示模拟结果，首先实现一维随机行走模型的可视化，每一次模拟中都包含两个醉汉的随机行走轨迹（由于是随机过程，所以每一次的模拟的结果都不相同，从这里我们就可以发现随机过程与混沌过程的差别所在）：
模拟一与模拟二的程序完全相同（并且每次模拟中两个醉汉的执行的是同一个程序），结果却差别很大。而对于混沌过程而言，只要能够确保初始条件完全相同，是能够保证两次程序运行的结果完全相同的，只是结果对于初值条件的依赖性极其敏感而已。
虽然我的计算机还能够继续模拟下去，但是一维情况的结果并不令人惊艳，所以一维随机行走的可视化到此为止。最后插入一个色彩绚丽令人眼前一亮的折线图表示法：
（P.S. 来自《维基百科》，但是窃以为程序十分简单，每个人都可以画出来）
下一步，我想要研究的是醉汉系综行走n步之后离原点的位移与n的关系，我分别采用了500和5000个醉汉所组成的系综取平均值并对结果进行了线性拟合，以此来观察醉汉<x>与n的线性关系（每一幅图的下方都给出了拟合的函数结果）：
<x> = 0.0005352 n + 0.03009
<br >
<x> = - 0.0001716 n - 0.02086
所以从上面两幅图的模拟结果中，我们很容易得出结论，即对于一维随机行走过程而言，<x>是趋向于零的，并且从图中点的在x=0附近的分布趋势可以看出，n越大，点的分布越弥散，代表方差越大；而系综越大（醉汉越多），则点的分布越集中，代表方差越小，这一点可以由概率论与数理统计的中心极限定理予以证明，在此就不赘述了。
再下一步，我所研究的就是醉汉系综行走了n步之后到原点距离的平方与n的关系，我分别采用了500，5000和50000个醉汉所组成的系综取平均值并对结果进行了线性拟合，以此来观察醉汉<x²>与n的线性关系（每一符图下方都给出了拟合的函数结果）：
<x²> = 0.9838 n - 0.8176
<x²> = 0.9952 n - 0.7888
<x²> = 1.000 n - 0.9740
从图中我们可以得出与上面相同的结论，不过这一次我们所验证的公式为<x²> = n（即扩散定律中<x²> = 2Dt中D取1/2）。
定步长，不等概率
下面我们考虑左右不等可能随机行走，即醉汉具有一定的倾向性，有更大的概率走向左边或者是右边，我们不妨考虑其中的一种情况即可，假设醉汉每一步向左走单位长度的概率是3/4，向右走单位长度的概率是1/4。那么，我们只需要将产生的随机数标尺的比较值从1/2改为3/4即可。模拟的结果为（依照惯例，我分别给出了100，1000和10000步的模拟结果）：
从图中我们可以看出，在具有倾向性的情形中，两个醉汉系综的轨迹相似性提高了很多，这也意味着，两个醉汉系综的位移之差的平均值减小了许多，并且，我们也可以预料到，倾向性越明显，轨迹的相似性就越高，位移之差的平均值就越小。这些结论，运用最基本的概率论知识都可以轻松地证明（求出解析表达式）。
不定步长，等概率
假设醉汉每一步所行走的长度为是不确定的，但是一定属于0到1之间，那么，我们可以用0到1之间的随机数来模拟醉汉每一步的步长，或者直接省去随机数标尺的步骤，将醉汉的步长设置为-1到1之间的随机数，似乎更为简洁。然后，我们重复之前的步骤，首先是可视化的轨迹结果（按照惯例，我分别给出了100，1000和10000步的模拟结果）：
第三幅图10000步我只给出了一个系综的模拟结果，这是为了图片的清晰起见。从上面的三幅图中，我们可以看出，在不定步长的情况下，醉汉的行动范围较之定步长的情况发生了明显的缩小，这也是符合我们预期的，不过，我们真正关心的是<x²>与n的关系，对其的模拟结果如下（按照惯例，我们分别采用500和5000个醉汉所组成的系综，并在每一符图的下方我们照旧给出拟合的函数结果）：
<x²> = 0.3340 n - 0.3906
<x²> = 0.3329 n - 0.3226
此时，与定步长的情形不同，<x²>不再与n相等，它们的关系依旧是线性关系，只不过系数发生了变化，新的关系式为<x²> = n/3, 所对应的扩散定律中的D为1/6。
接下来我要进行的是二维与三维随机行走的模拟，所研究的问题与一维相同，甚至结论也大体相同，所以，之后的我的文字叙述将更加简洁，只强调与一维情况有所不同的地方。
二维随机行走的模拟
二维随机行走不同于一维的情况，我把它分为两种情况来研究，一种是lattice random walk（格点随机行走），顾名思义，即醉汉每一步只有四个可选方向：上，下，左，右。另一种情况是全方向随机行走，顾名思义，即醉汉每一步都可以选择平面内任意一个方向。对于二维随机行走，我将只讨论等概率的情况。
定步长，lattice
二维定步长lattice random walk与一维左右等可能随机行走唯一的区别在于增加了两个可选的方向，所以我需要做出的改变仅仅在于将随机数标尺分别于1/4,1/2,3/4进行比较，然后选择醉汉下一步行走的方向。按照惯例我们首先给出可视化的轨迹模拟，下面是一个醉汉二维定步长lattice random walk随机行走3步的所有可能情况中的三种：
接下来分别是是两个醉汉二维定步长lattice random walk随机行走100，1000，10000,100000和1000000步的轨迹（用黑色实心圆标注出了每个醉汉的起点和终点）：
然后探究<ρ²>与n的关系，按照惯例，分别采用500和5000个醉汉所组成的系综（照旧，图下给出拟合所得的函数结果）：
<ρ²> = 1.003 n - 1.205
<ρ²> = 0.9944 n - 0.7177
所得的结果与一维的完全相同。
定步长，全方向
对于二维定步长全方向随机行走，我们可以采用极坐标的参数来实现全方向均匀分布，对于极坐标而言即是θ在0到2π之间均匀分布，ρ = 1。所以我将θ取为0到2π之间的随机数，然后利用极坐标与直角坐标的关系x = ρcosθ， y = ρsinθ，将极坐标转化为直角坐标再进行作图（也可以选择直接在极坐标系中作图，Matplotlib可以实现此功能），所得结果如下，按照惯例，首先给出可视化的轨迹模拟图，以下分别是两个醉汉二维定步长全方向随机行走100,1000,10000，100000,1000000步的结果：
我注意到，当步数足够多时，仅从所给出的轨迹图，我们几乎已经无法区分lattice random walk与全方位随机行走（都是二维定步长）了。当然，将如果可以将图的比例尺放大，它们仍然是十分容易区分的。还有一种区分的方法，即为观察醉汉行走的范围，很明显lattice random walk大于全方位随机行走。
另外，我给出了另一种画法（在折线标示出醉汉轨迹的基础上用黑点标注出醉汉在第n步的位置），以下的结果分别是一个醉汉二维定步长全方位随机行走100,1000和10000步的轨迹图：
若以白色的折线标注醉汉的轨迹，10000步的情况如下图：
接下来仍然是它的<ρ²>与n的关系：
<ρ²> = 1.003 n - 1.205
<ρ²> = 0.9944 n - 0.7177
结果与二维定步长lattice random walk相同，也与一维定步长的情况相同。
不定步长，全方向
二维不定步长全方向随机行走意味着在二维全方位随机行走的基础上，醉汉每一步的步长都是不确定的，但是其大小应为0到1之间的均匀分布。所以我们可以采用极坐标，将ρ取为0到1的随机数；或者直接采用直角坐标，将醉汉每一步的位移正交分解为x，y两个方向，在每个方向的位移取为从负二分之根号二到二分之根号二的随机数。那么，我们的结果是，首先，一个醉汉，三步，三种情况的轨迹图：
下面是，两个醉汉，100,1000,10000,100000和1000000步的轨迹图：
当行走的步数足够多之后，我们也只能通过醉汉行走的范围来区分定步长lattice random walk，全方向随机行走与不定步长的随机行走，他们行走范围的大小关系排序为（一般）：定步长lattice random walk > 定步长全方位随机行走 > 不定步长随机行走。
然后是<ρ²> 与n的关系，如下:
<ρ²> = 0.3324 n - 0.3724
<ρ²> = 0.3332 n - 0.3190
结果与一维不定步长随机行走相同。
接下来我会介绍三维随机行走的模拟，由于其<r²>与n的关系与在一维、二维随机行走相对应的情况下完全相同，所以省略。
三维随机行走的模拟
三维随机行走是最有趣的一部分，三维随机行走也可以分为类似二维的三维lattice randon walk与三维全方位随机行走但是如何实现三维全方位随机行走却费了作者的一番周折。在探究的过程中，作者也犯过一个典型的错误，即直接采用球坐标的两个参数来操作，我将错误的取法（空间分布并不均匀）得出的结果也列出，与正确的结果进行比较，大家就通过观察图中差异发现错误的取法错在哪里。
定步长，lattice random walk
按照惯例：
定步长，全方位随机行走 - 正确做法
我为大家简单地介绍一下三维球面上的Marsaglia方法，这是一种基于抽样变换的方法，产生的结果非常漂亮。原理在此我就不讲了，只给出实行的方法。
第一步：随机抽样产生一对-1到1之间均匀分布的随机数u,v。
第二步：计算r² = u² + v², 如果r² > 1，则重新抽样，直到满足r² < 1。
第三步：计算
即得到球面上均匀分布的点，效果图是这样的，360°全方向均匀无死角：
以上两张图片摘自博客园
按照惯例，模拟结果为：
我们会发现醉汉在x,y,z三个维度的行走范围是大概相当的，这也证实了我们采用三维球面上的Marsaglia方法产生全方向随机行走的正确性。
定步长，全方位随机行走 - 典型错误
一种典型的错误做法是套用二维极坐标的思路，在三维中使用球坐标的参数直接产生随机角度，从而生成随机方向分布，具体做法为：
第一步：分别生成-π/2到π/2的随机数θ和0到2π的随机数φ。
第二步：通过球坐标与直角坐标的关系，计算
但是，这样产生的其实并不是随机分布，究其原因，这是因为立体角微元的表达式|dΩ| = |sinθdθdφ|,不是与θ和φ的线性关系。所以，从这个表达式中，我们可以看出，当sinθ越大时，即越接近于±π/2时，在该方向取点的可能性越大；反之，若sinθ越小时，即越接近于0时，在该方向取点的可能行越小，反映到图像上的结果也十分明显，我们可以看到球面上越靠近两极处点的密度越大，越靠近赤道处点的密度越小，结果如下：
以上两张图片同样摘自博客园
按照惯例，模拟结果为：
从图中我们可以看出，很明显的，醉汉在z方向上的行走范围与x，y方向上的行走范围根本不是一个数量级的，在z方向上的行走范围远远大于在x，y方向上的行走范围，由此我们也可以从侧面看出这种做法并不能实现全方向随机行走，而是具有明显的两极倾向性。
此时<r²>与n的关系发生了变化，不再是线性的了。于是我采用了二次多项式进行拟合，并且在图下给出了拟合结果的函数关系。
<r²> = 0.4069 n² - 0.3572 n + 1.310
<r²> = 0.4050 n² - 0.1890 n + 0.5491
不定步长，全方向
按照惯例：
那么，我的图片展示环节到这里就结束了。
Ⅲ 结论（致谢）
具体的结论，我在正文的讨论中都已经提及，总结一下，整篇文章只不过是
实现了一维，二维和三维随机行走的轨迹模拟，其中，每一种维度都分为定步长与不定步长两类，此外，在二维与三维的情形中还分类讨论了lattice random walk与全方向随机行走两种随机行走方式。
并且分别验证了它们满足扩散定律,并且有着相同的扩散常数，唯一不同的是定步长与不定步长情形的扩散常数。
着重讨论了三维情形的全方位随机行走，简单地介绍了实现三维全方位随机行走的三维球面上的Marsaglia方法，其实，还有许多种方可以实现球面上的均匀分布，比如直接抽样法，力学方法（十分神奇，这种方法是基于力学原理产生的，球面上的点是主动运动的，他们相互排斥，直到一个最均匀的状态。思路是： 对球面上的每一个点施加 坐标： x,y,z; 施加速度，vx,vy,vz ; 计算其余点对当前点的作用力，求解一个微分方程组。初始状态随机分布，每一个计算步更新当前粒子的运动速度和位移。——摘自博客园），并且讨论了为什么生成两个随机分布的角度之后按照球坐标的形式产生的方法式不可行的。
这篇文章的所有内容均由作者独立完成。
程序代码
- 期末代码1
- 期末代码2
- 期末代码3
- 期末代码4
- 期末代码5
- 期末代码6
- 期末代码7
- 期末代码8
- 期末代码9
- 期末代码10
- 期末代码11
- 期末代码12
- 期末代码13
- 期末代码14
- 期末代码15
- 期末代码16
- 期末代码17
- 期末代码18
- 期末代码19
- 期末代码20
- 期末代码21
- 期末代码22
- 期末代码23
- 期末代码24
Ⅳ 论文引用
Giodano, N.J., Nakanishi, H. Computational Physics. Tsinghua University Press, December 2007
随机游走的无规则行走与扩散定律．中国百科网[引用日期2016-12-24]
欧阳钟灿. 攀登21世纪生命科学的通天塔——漫谈《生物物理学:能量、信息、生命》[J]. 科学（上海）, 2007, 第2期
Pearson, K. (1905). "The Problem of the Random Walk". Nature. 72 (1865): 294. doi:10.1038/072294b0.
Kohls (2016) Expected Coverage of Random Walk Mobility Algorithm, arxiv.org
齐民友等，《概率论与数理统计（第二版）》，高等教育出版社，August 2011.
