面向计算思维的蒙特卡罗C语言程序设计案例探究

顾丽红，丁淑妍

（中国石油大学（华东） 计算机与通信工程学院，山东 青岛 266580）

0 引 言

计算思维[1]概念由Jeannette M. Wing教授于2006年首次完整提出，随后引起了国内计算机教育界的重视。教育部计算机专业教学指导委员会明确提出将计算思维列入计算机专业人才的专业基本能力[2]。如何培养非计算机专业大学生的计算思维能力，这是摆在大学计算机基础教育工作者面前的难题。目前国内高校的C语言程序设计课程普遍存在注重知识基础而忽视能力基础的问题，能力的培养则需要引入更多能引起学生兴趣的程序设计案例实践[3-5]。近几年相关的教研论文越来越多，大多以数值计算类为主，面向跨计算学科的广义计算思维[6]的案例设计实践类论文还不多见。计算机模拟通常综合了计算机科学家和领域科学家的思维方式，在课程教学实践中引入这类程序设计案例，有利于培养学生的计算思维能力。

1 计算思维与计算机模拟

进入信息和数据时代，计算思维是每个人的基本技能，就像逻辑思维、数学思维、设计思维一样，是人们认识问题和解决问题的重要思维方式之一。计算思维能力是人们有意识地使用计算机科学家的思想、方法、技术、工具、资源、环境进行思考和活动的能力。如何全方位培养这种能力是教育工作者面临的挑战。

计算机模拟通常是指用抽象模型来模拟特定系统的计算机程序。人们在研究、设计、构造复杂系统时，往往需要设计制造一个模型来进行各种试验。传统方法是建立数学模型或实物模型，前者缺乏直观性也不便于试验，后者比较直观但不够经济和方便。随着计算机的出现和处理能力的不断增强，利用计算机来模拟各种复杂系统，成为重要的科学研究手段。计算机模拟的一般过程如图1所示。

其中，①数学建模与形式化。明确模拟工作的目标，确定模拟方案的评价准则。②模拟建模。根据问题特点和模拟要求从稳定性、精度和性能3个维度设计合适的算法。③程序设计。把模拟模型用计算机能执行的程序来描述。④模拟运行。分析运行结果是否合适，通过前几步查找问题，不断修正，直到结果满意。⑤模拟实验和结果处理。每一个方案，采用不同随机数序列重复多次，采用数理统计方法分析模拟结果的统计特征；不同方案，采用相同随机数序列进行模拟运行，消除因随机数序列不同而引起的差异；模拟结果图形化、可视化，乃至虚拟现实环境。

图1 计算机模拟的一般处理过程

计算机模拟的设计原则包括：①分级模拟原则。采用与任务匹配的模型或算法，合理简化问题，突出问题的关键。②效率兼顾原则。在模拟的不同阶段，合理调整对准确度和速度的要求，提高模拟效率。③可信验证原则。所有模拟结果应该有可信的验证方法或依据。

在教学实践中，让学生充分理解计算机模拟过程中的算法构建、迭代修正、数理统计方法，以及结果图形化和可视化的重要性，在提升学生兴趣的实践中培养能力。

2 蒙特卡罗方法思想

蒙特卡罗方法是一种随机模拟方法，以概率和统计理论方法为基础，使用随机数（伪随机数）解决计算问题，又称为统计模拟法。因其具有概率统计特征而获得数学家冯·诺依曼用赌城蒙特卡罗来命名。蒙特卡罗方法起源于1777年法国数学家蒲丰（Buffon）提出用投针试验计算圆周率π，进而推广到用概率法计算不规则图形面积，甚至用小规模抽样调查进行民意测验，来预测竞选的优胜者，也是同样的思想。当然科学计算中的问题比这复杂很多，比如金融领域的交易风险评估，问题的维数（随机变量的个数）高达数百、几千，问题的难度成指数增长，传统数值方法难以胜任，而蒙特卡罗方法因计算复杂性不依赖于维数则可以很好地对付维数灾难。

蒙特卡罗方法具有简单、快速、节省存储单元的优点，省却了一般人难以掌握的复杂的数学推导和演算过程，而且不受问题的几何形状复杂性影响，适用于处理大型复杂问题。随着计算机处理复杂系统的能力日益增强，蒙特卡罗方法的应用也越来越广泛。它不仅较好地解决了多重积分计算、微分方程求解、积分方程求解、特征值计算和非线性方程组求解等高难度和复杂的数学计算问题，而且在计算物理学、核物理、金融工程学、宏观经济学、生物医学、可靠性及计算机科学等广泛的领域都得到成功的应用。

在教学实践中，通过引入蒲丰投针问题和不规则面积计算问题，学生对蒙特卡罗方法思想有了浓厚的兴趣，为之后的C语言算法设计思想的养成打下基础。

3 蒙特卡罗模拟案例设计

蒙特卡罗模拟的基本思想：当所求问题的解是某个事件的概率，或是某个随机变量的数学期望，或是与概率、数学期望有关的量时，通过某种试验的方法，得出该事件发生的频率或者该随机变量的若干个具体观察值的算术平均值，进而得到问题的解。采用蒙特卡罗方法进行计算机模拟的步骤：首先设计反映系统各部分运行时逻辑关系的逻辑框图（模拟模型），然后通过具有各种概率分布的模拟随机数来模拟随机现象。

3.1 3种蒙特卡罗π的计算模拟

法国数学家蒲丰于18世纪提出投针问题：设有一个以平行且等距木纹铺成的地板，现在随意抛一支长度比木纹之间距离小的针，求针和其中一条木纹相交的概率，并以此概率可以近似计算圆周率。

蒲丰投针问题的数学建模：平行线距离为2a，针的长度为2l，且l≤a。设针投到地面上是随机的，所以位置可以用二维随机变量(x,θ)来描述，x为针中心的坐标，θ为针与平行线的夹角。任意投针，就是意味着x与θ都是任意取的，但x的范围限于[0,a]，夹角θ的范围限于[0, π]。在此情况下，针与平行线相交的数学条件是：x≤l×sinθ。x和θ均服从均匀分布，且相互独立，通过(x,θ)的概率密度函数求出概率分布：P= 2l/πa，所以π的近似值为：2l/aP。

蒲丰投针C语言算法设计思路：

（1）随机生成一个介于[0,a]表示针中心坐标的数x和一个介于[0, π]表示针与平行线的角度y；

（2）如果x≤l×siny，则表示针与平行线相交；

（3）计算π的近似值。

C语言程序设计实现如下：

蒲丰投针实验的重要性在于它是第一个用几何形式表达概率问题的例子，开创了使用随机数处理确定性数学问题的先河。在教学实践中，蒲丰实验问题新颖奇妙，引起学生极大的兴趣。为了达到举一反三的效果，案例设计中引出蒙特卡罗思想计算π的另外两种解法：圆中投点（几何法）和随机数互质法。

圆中投点问题描述：有一正方形木板，以及内切正方形内接圆盘，随机向正方形中投点，求点落在圆盘中的概率，并以此概率近似计算圆周率。

圆中投点数学建模：有一个以（0, 0）为中心的边长为2的正方形，以及这个正方形中半径为1的内接圆，随机向正方形中投点，求点落在内接圆中的概率。设点投到木板上是随机的，所以位置可以用二维随机变量(x,y)来描述。任意投点，就是意味着x与y都是任意取的，二者的范围限于[0, 1]。在此情况下，落在内接圆中的数学条件是：x2+y2≤ 1。x和y均服从均匀分布，且相互独立，求出概率分布：P= π/4，所以π的近似值为：4×P。

圆中投点C语言算法设计思路：

（1）随机生成两个介于[0, 1]区间的表示点坐标的x和y；

（2）如果x2+y2≤1，则表示点落在内接圆中；

（3）计算π的近似值。

C语言程序设计实现如下：

随机数互质问题描述：R·查特在1904年发现，两个随意写出的数中，互质的概率为6/π2，并以此概率近似计算圆周率。

随机数互质数学建模：任意产生两个整数x和y，x和y均服从均匀分布，且相互独立，通过数学推导（推导过程较复杂）可以求出x和y互质的概率为：P= 6/π2，所以π的近似值为6/P的算术平方根。

随机数互质C语言算法设计思路：

（1）随机生成两个介于[0, 32 767)的x和y；

（2）如果x和y互为质数，则统计变量累加1（注意特别处理生成数0的情况）；

（3）计算π的近似值。

C语言程序设计实现如下：

测试实验的硬件环境：Intel i5 M540处理器，主频2.53GHz，双核四线程，3MB高速缓存。软件环境：Window7旗舰版，CodeBlock 13.12 IDE环境，TDM-GCC v4.8.1编译器。为了尽可能减少其他应用程序对测试的干扰，测试前尽可能关闭开机启动项的应用软件，比如杀毒软件、下载工具等，并禁用网络。采用多次测试取均值的方法可以降低测试误差，本实验取10次测试的平均值，实验结果如表1。

众所周知，圆周率π的真值在3.141 592 6和3.141 592 7之间，假设以中间值3.141 592 65为真值，我们发现当模拟的次数N为千万次和亿次的误差分别为： 0.438‰和0.048‰， 0.235‰和0.023‰， 0.071‰和0.030‰，π的精度有了明显的提高；当N为10亿次时，精度没有明显提升，求更高精度的π需要寻求其他方法。蒙特卡罗方法的收敛速度与N的算术平方根成比例，意味着要使结果精度提高一位，应该增加一百倍的模拟计算工作量。

3.2 4种概率分布的随机变量模拟

蒙特卡罗方法需要模拟随机事件，即需要随机数，那么计算机是如何产生随机数的呢？随机数是来自统计学的概念，真正的随机数（真随机数）是使用物理现象产生的，比如常见的彩票的摇号、掷钱币、骰子、转轮等，密码学领域需要真随机数。目前计算机产生的随机数，是通过一个固定的、可以重复的算法产生，不是真正的随机，但是具有类似于随机数的统计特征，通常称为伪随机数，可以满足蒙特卡罗方法的要求，所以默认计算机产生的随机数是随机的。

用蒙特卡罗方法模拟系统或过程时，需要不同概率分布的随机数，通常先生成均匀分布的（0, 1）随机数，这是其他分布随机数的基础。理论上说，具有连续分布的随机数，通过数学变换或近似等方法，可以生成其他任意分布的随机数。设μi(i=1, 2, …)是区间[0, 1]内的均匀分布的独立随机变量，而另一给定分布函数F(x)的随机变量为则这一随机变量xi可以由其反分布函数求得其抽样值，即

均匀分布随机变量算法设计思路：已知连续型随机变量在有限区间内取值，则其概率密度函数为其分布函数为知所以

指数分布随机变量算法设计思路：已知连续型随机变量在有限区间内取值，则其概率密度函数为故其分布函数为知若(1-μ)是（0, 1）均匀分布随机数，则可用下式简化为

正态分布随机变量算法设计思路：已知连续型随机变量服从正态分布N(μ, σ2)，其概率密度函数为取区间（0, 1）上两个均匀分布随机数μ1和μ2，利用二元函数变换可得到个独立的标准正态分布随机变量。

表1 三种蒙特卡罗π的计算模拟结果

泊松分布随机变量算法设计思路：泊松分布源自二项分布，在二项分布的伯努利试验中，如果试验次数n很大，二项分布的概率p很小，且乘积λ=np比较适中，则事件出现的次数的概率可以用泊松分布来逼近。泊松分布是离散型概率分布，表示固定尺度的连续区间上给定的事件发生次数的概率，其中n可以看作无穷大，泊松分布单位时间内随机事件发生的次数满足泊松分布，比如电话交换机接收到呼叫的次数、汽车站台的候客人数、机器出现故障的次数、自然灾害发生的次数、DNA序列的变异数等。

在C语言程序设计中，利用计算机产生的均匀分布随机数，通过概率积分变换算法可以比较容易得到满足其他概率分布的随机变量。C语言中产生随机数，通常用到两个函数：srand()和rand()。srand()用来为计算机产生随机数设置seed种子，否则每次程序运行产生的随机数序列是一样的。避免这种情况的通常方法是采用计算机的当前时间作为seed种子，因为时间是在不断变化的。rand()函数随机生成满足均匀分布的[0,RAND_MAX]（RAND_MAX 为 0x7fff，即 32 767）之间的整数。如果要使范围大一点，可以通过产生几个随机数的线性组合来实现任意范围内的均匀分布随机数，而且经过一定的四则运算和取模运算，可以比较容易地得到任意区间的随机变量，比如：均匀分布、指数分布、正态分布和泊松分布。

产生高精度均匀分布随机数。用rand()函数产生随机数，需要对精度有所认知。rand()的随机数分辨率为32 767，两个也就是65 534，如果需要要产生（-1 000, 1 000）之间且精度为4位小数点的均匀分布，则分辨率要求为1 000×10 000×2=20 000 000，这样显然远远不够，但可以用两个随机数的乘法来达到精度要求，C语言程序设计片段如下。

利用[0, 1]区间的均匀分布可以产生指数分布随机数，C语言程序设计片段如下。

利用[0, 1]区间的均匀分布可以产生正态分布随机数（Box Muller方法），C语言程序设计片段如下。

利用[0, 1]区间的均匀分布可以产生泊松分布随机数，Knuth算法设计思路：

1）l赋初值e-λ，k赋初值0，p赋初值1.0；

2）生成(0, 1)区间均匀分布随机数u，p赋值为p*u，k累加1；

3）如果p>l，则重复2），直至不满足条件；

4）返回k-1。

C语言程序设计片段如下：

实验软硬件环境同上，实验结果如图2（频率曲线图）。

从4种分布的频率曲线图可以看出，通过线性组合和四则运算，计算机产生的伪随机数可以满足蒙特卡罗方法的要求，在实际生产、生活和工程项目中可以大量应用。

4 结 语

蒙特卡罗模拟有很多专门的工具软件，比如Mathmetica、MatLab、SPSS、CrystalBall等，其功能强大且容易使用。本文仅从C程序设计课程教学和实践环节的需要，提出通过计算机模拟随机事件引导学生对计算思维的认知和培养学生动手实践的能力。大数据时代需要培养更多理解概率论和数理统计原理的融合型复合人才，建议大学计算机基础课程教学工作者根据学生的专业背景设计更多跨学科领域且体现计算思维思想的案例库进行交流和分享。下一步将在本文工作的基础上，从大数据分析的角度，研究设计诸如投资风险分析、随机库存预测等进一步培养计算思维理能力的综合型教学案例。

图2 计算机生成的4种分布随机数的频率曲线图

[1]Wing J M. Computational thinking[J]. Communications of the ACM, 2006, 49(3): 33-35.

[2]教育部高等学校计算机科学与技术教学指导委员会. 高等学校计算机科学与技术专业人才专业能力构成与培养[M]. 北京: 机械工业出版社, 2010.

[3]顾丽红, 李传秀, 吴少刚. 培养计算思维能力的矩阵乘法C语言程序设计案例探究[J]. 计算机教育, 2016(1): 149-152.

[4]姚天昉. 在程序设计课程中引入“计算思维”的实践[J]. 中国大学教学, 2012(2): 61-62.

[5]陈文智, 陈越, 庄越挺. 面向系统设计能力培养的教学改革探索[J]. 计算机教育, 2013(20): 70-76.

[6]蒋宗礼. 计算思维之我见[J]. 中国大学教学, 2013(9): 5-10.