1、10:22:54MCM1第九章、随机模型第九章、随机模型 9.3 9.3 拉丁方与正交设计拉丁方与正交设计 9.1 9.1 几何概率模型几何概率模型 9.2 9.2 计算机模拟计算机模拟10:22:54MCM2 现实世界中充满了不确定性，我们所研究的对象往往现实世界中充满了不确定性，我们所研究的对象往往受到诸多随机因素的影响，因此所建立的数学模型涉及到受到诸多随机因素的影响，因此所建立的数学模型涉及到的变量往往并非确定型变量而是随机变量，甚至有时候变的变量往往并非确定型变量而是随机变量，甚至有时候变量之间的关系也不是确定的函数关系，我们将这类模型称量之间的关系也不是确定的函数关系，我们将这类模
2、型称为随机模型。为随机模型。9.1 9.1 几何概率模型几何概率模型 一类涉及到一类涉及到“等可能性等可能性”的概率问题我们称之为几何的概率问题我们称之为几何概型或几何概率模型，我们先来观察一下下面的例子。概型或几何概率模型，我们先来观察一下下面的例子。10:22:54MCM3例例9.1 9.1 约会问题约会问题 一位男生与一位女生约定晚饭后一位男生与一位女生约定晚饭后1818时到时到1919时之间在教时之间在教学大楼后的左边第三棵柳树下会面。双方约定，先到者必学大楼后的左边第三棵柳树下会面。双方约定，先到者必须等候另一个人须等候另一个人1515分钟，过时如另一人仍未到达则可离去，分钟，过时如
3、另一人仍未到达则可离去，问两人会面的机会有多大？问两人会面的机会有多大？显然，男生和女生到达约会地点的时间是随机的、不显然，男生和女生到达约会地点的时间是随机的、不确定的，故所涉模型应为随机模型。根据题意，我们已经确定的，故所涉模型应为随机模型。根据题意，我们已经知道了男生和女生各自到达的时间范围，他们均可能在知道了男生和女生各自到达的时间范围，他们均可能在1818时到时到1919时之间的任何时刻到达约会地点。但是无法预知他时之间的任何时刻到达约会地点。但是无法预知他们到达的确切时刻，那么等候另一人的们到达的确切时刻，那么等候另一人的1515分钟自然也不确分钟自然也不确定，问题好像有点麻烦。定
4、，问题好像有点麻烦。10:22:54MCM4 如果到达的时间只有有限多个时间点，那么所涉及的如果到达的时间只有有限多个时间点，那么所涉及的数学模型属于古典概型，我们只需列出男生和女生到达时数学模型属于古典概型，我们只需列出男生和女生到达时间的所有可能组合，然后从中找出所有可以会面的组合，间的所有可能组合，然后从中找出所有可以会面的组合，用可以会面的组合数去除以所有可能的组合数就可得到两用可以会面的组合数去除以所有可能的组合数就可得到两人会面的概率。但现在到达时间可取自一个无限集合人会面的概率。但现在到达时间可取自一个无限集合(区区间间)，这给我们的计算带来了一些困难。用古典概型的方，这给我们的
5、计算带来了一些困难。用古典概型的方法无法解决此问题，我们要另寻其它方法。法无法解决此问题，我们要另寻其它方法。几何概型讨论的是无限样本空间中的概率问题，在此几何概型讨论的是无限样本空间中的概率问题，在此空间中随机试验的每一结果都是等可能发生的。在几何概空间中随机试验的每一结果都是等可能发生的。在几何概型中，随机试验的所有可能结果构成一个样本空间，样本型中，随机试验的所有可能结果构成一个样本空间，样本空间通常是一个几何区域空间通常是一个几何区域 。10:22:54MCM5 试验中可能发生的事件则为试验中可能发生的事件则为 中中的一个子区域，而随机事件发生的概的一个子区域，而随机事件发生的概率则是
6、在对两者进行了比较之后计算率则是在对两者进行了比较之后计算出来的。比如我们在一个面积为出来的。比如我们在一个面积为 的的区域区域 中，等可能地任意投点中，等可能地任意投点(如图如图9.1)9.1)。在区域。在区域 中有一个小区域，它中有一个小区域，它的面积为的面积为 ，投点落入，投点落入A A中的可能性中的可能性大小只与大小只与A A的面积的面积 成正比，成正比，图图9 91 1而与的位置及形状无关。如果仍将而与的位置及形状无关。如果仍将“投点落入小区域投点落入小区域A A”这个随机事件记作这个随机事件记作A A，并设，并设 (即投点必落在中即投点必落在中)，可，可得事件得事件A A发生的概率
7、为发生的概率为 ，这样的概率问题通常，这样的概率问题通常被称为被称为”几何概型几何概型”。SASAS()1P()ASP AS10:22:54MCM6 应当注意的是，这里的面积只是二维平面内某个区域应当注意的是，这里的面积只是二维平面内某个区域的测度而已，如果是在一条线段上等可能地任意投点，公的测度而已，如果是在一条线段上等可能地任意投点，公式里的面积应改为长度；如果在一个三维立方体内等可能式里的面积应改为长度；如果在一个三维立方体内等可能地投点，则面积应该为体积。地投点，则面积应该为体积。现在，我们来求解例现在，我们来求解例9.19.1。设。设x x和和y y分别表示男生和女分别表示男生和女生
8、到达约会地点的时间生到达约会地点的时间(为计算方便，从为计算方便，从1818时开始计时时开始计时)，建立平面直角坐标系如图建立平面直角坐标系如图9-29-2所示，所有可能到达时间的所示，所有可能到达时间的组合，即组合，即(x,y)(x,y)的所有可能结果构成一个边长为的所有可能结果构成一个边长为6060的正方的正方形形(以分钟为单位以分钟为单位)。另外，由题意两人能够会面的充要条。另外，由题意两人能够会面的充要条件是件是15xy10:22:54MCM7图图9-29-2 01020304050600102030405060 xyA 10:22:54MCM8 可能会面的时间组合由图中的阴影部分所表
9、示。我们可能会面的时间组合由图中的阴影部分所表示。我们假设两人到达约会地点的时间在这一小时中均是等可能的，假设两人到达约会地点的时间在这一小时中均是等可能的，此时，例此时，例9.19.1就成为了一个几何概型问题。记二人会面的就成为了一个几何概型问题。记二人会面的事件为事件为A A，由等可能性可知二人会面的可能性为，由等可能性可知二人会面的可能性为22260457()6016ASP AS 几何概型的用途十分广泛，其中最为著名的是利用它几何概型的用途十分广泛，其中最为著名的是利用它来求来求 的近似值。在的近似值。在17771777年出版的年出版的或然性算术实验或然性算术实验一一书中，法国科学家蒲丰
10、（书中，法国科学家蒲丰（BuffonBuffon）提出了著名的蒲丰投针）提出了著名的蒲丰投针试验，并以该实验方法计算试验，并以该实验方法计算 。10:22:54MCM9 这个实验方法的操作很简单：在一张白纸上画出一组这个实验方法的操作很简单：在一张白纸上画出一组间距为间距为 的平行线，并找一根粗细均匀，长度为的平行线，并找一根粗细均匀，长度为 的细针，然后一次又一次地将小针任意投掷在白纸上。这的细针，然后一次又一次地将小针任意投掷在白纸上。这样反复地投掷多次，记录下针与这组平行线中的任意一条样反复地投掷多次，记录下针与这组平行线中的任意一条相交的次数，就可以得到相交的次数，就可以得到 的近似值
11、。例如，在某次实验的近似值。例如，在某次实验中，他选取了中，他选取了 ，然后，共计投了，然后，共计投了22122212次，其中针次，其中针与平行线之一相交了与平行线之一相交了704704次，这样求得圆周率的近似值为次，这样求得圆周率的近似值为 2212/704=3.1422212/704=3.142。下面，我们来简单地介绍一下该试验。下面，我们来简单地介绍一下该试验的原理。的原理。a()l la/2la10:22:54MCM10 如图如图9-39-3所示，用所示，用x x表表示针的中点与最近一条平示针的中点与最近一条平行线间的距离，用行线间的距离，用 表示表示针与此直线间的夹角，显针与此直线间
12、的夹角，显然有然有020ax图图9-39-310:22:54MCM11 由这两个不等式可以确定由这两个不等式可以确定 平面上的一个矩形区平面上的一个矩形区域域 (见图见图9-4)9-4)，其面积为：，其面积为：x()2aS 图图9-49-4为使针与平行线相交，其为使针与平行线相交，其充要条件是：充要条件是：sin20lx10:22:54MCM12这两个不等式确定的平面区域这两个不等式确定的平面区域A A的面积为的面积为由细针落纸的等可能性知，细针与平行线相交的概率由细针落纸的等可能性知，细针与平行线相交的概率 0()sin2lS Adl()2()S AlpSa 记记N N为投针试验中的试验次数
13、，为投针试验中的试验次数，n n为针与平行线相交为针与平行线相交的次数，则的次数，则n/Nn/N为细针与平行线相交的频率。由概率论中为细针与平行线相交的频率。由概率论中的大数定理，当试验次数足够多的大数定理，当试验次数足够多(N(N足够大足够大)时，可将事件时，可将事件发生的频率看作事件发生的概率的近似值，即发生的频率看作事件发生的概率的近似值，即于是得到于是得到 。2lnaN2lNan10:22:54MCM13 历史上有许多学者曾亲自做过这个试验，下表记录历史上有许多学者曾亲自做过这个试验，下表记录了他们的试验结果（把了他们的试验结果（把 折算成单位长度折算成单位长度1 1）:a 表表9.1
14、9.1试验者试验者年份年份 针长针长投掷次投掷次数数相交次数相交次数 得到得到 的近似值的近似值WolfWolf185018500.80.850005000253225323.15963.1596SmithSmith185518550.60.6320432041218.51218.53.15543.1554DeMorgan,C.DeMorgan,C.186018601.01.0600600382.5382.53.1373.137FoxFox188418840.750.75103010304894893.15953.1595LazzeriniLazzerini190119010.830.8334
15、083408180818083.14159293.1415929ReinaReina192519250.54190.5419252025208598593.17953.179510:22:54MCM14 蒙特卡罗蒙特卡罗(Monte Carlo)(Monte Carlo)方法，也称为计算机随机模方法，也称为计算机随机模拟方法，是一种基于拟方法，是一种基于“随机数随机数”的计算方法。这一方法源的计算方法。这一方法源于美国在第一次世界大战期间研制原子弹的于美国在第一次世界大战期间研制原子弹的“曼哈顿计曼哈顿计划划”。该计划的主持人之一，数学家冯诺伊曼用驰名世界。该计划的主持人之一，数学家冯诺伊曼用
16、驰名世界的赌城的赌城摩纳哥的摩纳哥的Monte CarloMonte Carlo来命名这种方法，使它来命名这种方法，使它蒙上了一层神秘的色彩。蒙上了一层神秘的色彩。9.2 9.2 计算机模拟计算机模拟 Monte Carlo Monte Carlo方法的基本思想源于我们上节介绍过的方法的基本思想源于我们上节介绍过的蒲丰投针试验。早在蒲丰投针试验。早在1717世纪，人们就已经知道用事件发生世纪，人们就已经知道用事件发生的的“频率频率”作为事件的作为事件的“概率概率”的近似值。只要设计一个的近似值。只要设计一个随机试验，使一个事件的概率于某未知数有关，然后通过随机试验，使一个事件的概率于某未知数有
17、关，然后通过重复试验，以频率近似表示概率，即可求得该未知数的近重复试验，以频率近似表示概率，即可求得该未知数的近似解。似解。10:22:54MCM15 显然，利用随机试验求近似解，试验次数要相当多才显然，利用随机试验求近似解，试验次数要相当多才行，比如前述的各位学者为了求，都亲自重复投针达几千行，比如前述的各位学者为了求，都亲自重复投针达几千次之多，这显然既费时又费力，非常麻烦。随着本世纪次之多，这显然既费时又费力，非常麻烦。随着本世纪4040年代电子计算机的出现，人们便开始利用计算机来模拟所年代电子计算机的出现，人们便开始利用计算机来模拟所设计的随机试验，使得这种方法得到迅速的发展和广泛的设
18、计的随机试验，使得这种方法得到迅速的发展和广泛的应用。特别是近年来高速电子计算机的出现，使得用数学应用。特别是近年来高速电子计算机的出现，使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能。方法在计算机上大量、快速地模拟这样的试验成为可能。若对已知分布的随机变量进行模拟，需要产生一系列服从若对已知分布的随机变量进行模拟，需要产生一系列服从这种分布的随机数。这种分布的随机数。其中最基本的模拟是区间上的均匀分布随机数的产生。其中最基本的模拟是区间上的均匀分布随机数的产生。产生均匀分布的随机数的常用方法很多，有平方取中法、产生均匀分布的随机数的常用方法很多，有平方取中法、线性同余法等，这里我
19、们不作专门的介绍，因为很多软件线性同余法等，这里我们不作专门的介绍，因为很多软件均提供了生成随机数的专用命令，如均提供了生成随机数的专用命令，如MATLABMATLAB提供了提供了randrand函函数来生成上的均匀分布随机数等。数来生成上的均匀分布随机数等。10:22:54MCM16 计算机产生的随机数都是按照某种确定的算法产生的，计算机产生的随机数都是按照某种确定的算法产生的，它遵循了一定的规律，一旦初始值确定，所有随机数也就它遵循了一定的规律，一旦初始值确定，所有随机数也就随之确定，这显然不满足真正随机数的要求，因此我们称随之确定，这显然不满足真正随机数的要求，因此我们称这种随机数为伪随
20、机数。但只要伪随机数能通过独立性检这种随机数为伪随机数。但只要伪随机数能通过独立性检验、分布均匀性检验、参数检验等一系列的统计检验，就验、分布均匀性检验、参数检验等一系列的统计检验，就可以把它们当作真正的随机数那样放心地使用。可以把它们当作真正的随机数那样放心地使用。上节中的约会问题我们也可以通过计算机模拟来解决。上节中的约会问题我们也可以通过计算机模拟来解决。用用x x和和y y分别记男生和女生到达约会地点的时间分别记男生和女生到达约会地点的时间(仍从仍从1818时开时开始计时，但以小时为单位计算始计时，但以小时为单位计算)，由于二人均可能在，由于二人均可能在1818时到时到1919时之间的
21、任何时刻到达约会地点，且在这一小时之内的任时之间的任何时刻到达约会地点，且在这一小时之内的任何时刻到达的可能性均相同，所以何时刻到达的可能性均相同，所以x x和和y y均服从区间均服从区间(0,1)(0,1)上上的均匀分布。的均匀分布。10:22:54MCM17 取一个充分大的正整数取一个充分大的正整数n n，重，重复复n n次，每次独立地从区间次，每次独立地从区间(0,1)(0,1)中随机地取一对数中随机地取一对数x x和和y y，分别检，分别检验验|x-y|15/60|x-y|15/60是否成立。设是否成立。设n n次次试验中不等式成立的共有试验中不等式成立的共有m m次，则次，则二人会面
22、的可能性可近似地看成二人会面的可能性可近似地看成为为m/nm/n。图。图9.59.5为某次模拟时取为某次模拟时取10001000个随机数对得到的散点图，个随机数对得到的散点图，并由此得到二人会面的可能性大并由此得到二人会面的可能性大约 为约 为 0.4 4 4 00.4 4 4 0，与 实 际 数 值，与 实 际 数 值7/16=0.43757/16=0.4375的误差仅为的误差仅为0.00650.0065，如果我们希望据此方法得到更为如果我们希望据此方法得到更为精确的结果，可令精确的结果，可令n n取更大的数。取更大的数。00.20.40.60.8100.10.20.30.40.50.60.
23、70.80.91图图9.59.510:22:55MCM18 蒲丰投针试验也可进行计算机模拟。但一般只能用来蒲丰投针试验也可进行计算机模拟。但一般只能用来演示而不能真正用来求演示而不能真正用来求的近似值，因为产生服从均匀分的近似值，因为产生服从均匀分布的布的 的随机数时，需要知道均匀分布的区间范围的随机数时，需要知道均匀分布的区间范围（），但既然范围已经知道，我们也就无需再），但既然范围已经知道，我们也就无需再求求了！那么能否利用计算机，用模拟的方法来求了！那么能否利用计算机，用模拟的方法来求呢？呢？根据上节的几何概率我们也可以自行设计求根据上节的几何概率我们也可以自行设计求的方法。的方法。0,
24、02或，比如我们可以采取下面的方法：如图比如我们可以采取下面的方法：如图9-69-6所示，取一个所示，取一个二维数组二维数组(x,y)(x,y)，取一个充分大的正整数，取一个充分大的正整数n n，重复，重复n n次，每次次，每次独立地从区间独立地从区间(0,1)(0,1)中随机地取一对数中随机地取一对数x x和和y y，并分别检验，并分别检验 是否成立。设是否成立。设n n次试验中不等式成立的共有次试验中不等式成立的共有m m次，则次，则 221xy4mn10:22:55MCM19 图图9-69-6 10:22:55MCM20 大部分其它分布的随机数的产生均基于均匀的分布大部分其它分布的随机数
25、的产生均基于均匀的分布随机数，主要采用的方法有逆变换法等。逆变换法基于随机数，主要采用的方法有逆变换法等。逆变换法基于定理定理9.19.1。记随机变量记随机变量X X的分布函数为的分布函数为F(x)F(x)，定义它的逆函数，定义它的逆函数 ,y,y的取值范围为的取值范围为0 0y y1 1。1()inf:()Fyx F xy定理定理9.19.1 设随机变量设随机变量Y Y服从区间服从区间(0,1)(0,1)上的均匀分布，上的均匀分布，则随机变量则随机变量 的分布函数为的分布函数为F(x)F(x)。1()XFY证明证明 由分布函数的定义及上述函数由分布函数的定义及上述函数 的定义，可以的定义，可
26、以得到随机变量得到随机变量X X的分布函数为的分布函数为 1F1()()()P XxP FYxP YF x10:22:55MCM21 因为因为Y Y服从区间服从区间(0,1)(0,1)上的均匀分布，显然当上的均匀分布，显然当0 0y y1 1时总有时总有 ，故，故 ，证毕。，证毕。()P Yyy()()P XxF x 由定理由定理9.19.1，要产生分布函数为，要产生分布函数为F(x)F(x)的随机数，只要的随机数，只要先产生区间先产生区间(0,1)(0,1)上的均匀分布的随机数上的均匀分布的随机数y y，然后计，然后计算算 ，即为分布函数为即为分布函数为F(x)F(x)的随机数。合成和筛的随
27、机数。合成和筛选等方法也可以由均匀分布随机数产生其它分布的随机选等方法也可以由均匀分布随机数产生其它分布的随机数。某些特殊分布如标准正态分布数。某些特殊分布如标准正态分布N(0,1)N(0,1)、指数分布的、指数分布的随机数的产生还有专门的定理随机数的产生还有专门的定理(定理定理9.29.2为标准正态随机为标准正态随机数的产生定理数的产生定理)。大部分数学软件均提供了专门的函数来。大部分数学软件均提供了专门的函数来产生常用分布随机数，给我们带来了很大的便利，其原产生常用分布随机数，给我们带来了很大的便利，其原理大体就是这样的。理大体就是这样的。1()FY1()FY10:22:55MCM22定理
28、定理9.29.2 设设 、是独立同分布的服从区间是独立同分布的服从区间(0,1)(0,1)上上的均匀分布的随机变量，令的均匀分布的随机变量，令则随机变量则随机变量 与与 相互独立，且均服从标准正态分布。相互独立，且均服从标准正态分布。1X2X1211212212(2ln)cos(2)(2ln)sin(2)YXXYXX 1Y2Y 利用此定理可以同时产生两个服从标准正态分布的随利用此定理可以同时产生两个服从标准正态分布的随机数，首先产生两个区间机数，首先产生两个区间(0,1)(0,1)上的均匀分布的随机上的均匀分布的随机数数 、，然后计算，然后计算出出 ，则，则 和和 均为标准正态随机数。均为标准
29、正态随机数。1x2x1y2y12112(2ln)cos(2)yxx 12212(2ln)sin(2)yxx 10:22:56MCM23 蒙特卡罗方法现已被成功地运用到诸如解微分方程、蒙特卡罗方法现已被成功地运用到诸如解微分方程、解积分方程、求特征值、矩阵转置等许多传统的数学问解积分方程、求特征值、矩阵转置等许多传统的数学问题中去，尤其是常被用于计算多重积分。题中去，尤其是常被用于计算多重积分。让我们来看一下较为简单的数值积分让我们来看一下较为简单的数值积分计算一元计算一元函数定积分函数定积分 的蒙特卡罗方法的主要思想。在区的蒙特卡罗方法的主要思想。在区间间(a,b)(a,b)上随机选取上随机选
30、取n n个点，用这个点，用这n n个点上的函数值的平均个点上的函数值的平均值来代替值来代替f(x)f(x)在区间在区间(a,b)(a,b)上的平均值。由此思想可以得上的平均值。由此思想可以得到计算一元函数定积分近似值的公式到计算一元函数定积分近似值的公式()baf x dx11()()()nbiaif x dxbaf xn(1,2,)ix in其中其中 是区间是区间(a,b)(a,b)上均匀分布的上均匀分布的n n个随机数。个随机数。这个思想显然很容易推广到多重积分上去。这个思想显然很容易推广到多重积分上去。10:22:56MCM24对于二重积分，有相应的公式对于二重积分，有相应的公式1212
31、121()()(,)(,)nbdiiacibadcf x xdx dxf xxn 其中是其中是 区间区间(a,b)(a,b)上均匀分布的上均匀分布的n n个随机个随机数，数，是区间是区间(c,d)(c,d)上均匀分布的上均匀分布的n n个随机个随机数。数。1(1,2,)ix in2(1,2,)ixin例例9.29.2 计算以计算以 为底，以为底，以 为顶的曲顶柱体的体积，图为顶的曲顶柱体的体积，图9.79.7为该为该柱体的顶曲面的图形。柱体的顶曲面的图形。(,)|23,14Dx yxy23(1)zyxy10:22:56MCM25图图9-79-710:22:56MCM26 由微积分知识，该曲顶柱
32、体的体积为由微积分知识，该曲顶柱体的体积为 ，因此，因此问题转化为求该二重积分。但被积函数较为复杂，直接问题转化为求该二重积分。但被积函数较为复杂，直接求解析解较为困难，我们可以使用蒙特卡罗方法求其近求解析解较为困难，我们可以使用蒙特卡罗方法求其近似的数值解。步骤如下：似的数值解。步骤如下：Dzd(,)iixy(1)(1)给定一个较大的自然数给定一个较大的自然数n n，生成，生成n n对随机数对随机数 ，和和 分别是区间分别是区间(2,3)(2,3)和和(1,4)(1,4)上的均上的均匀分布的随机数。匀分布的随机数。(注：生成区间注：生成区间(0,1)(0,1)上的均匀分布随上的均匀分布随机数
33、机数u u后，令后，令 ，则，则x x是区间是区间(a,b)(a,b)上均匀上均匀分布的随机数分布的随机数)(1,2,)iy inix()xaba u10:22:56MCM27(2)(2)将将 和和 代入被积函数表达式中计算出相应的代入被积函数表达式中计算出相应的 。ixiyiz(3)(3)计算计算 ，即可得到曲顶，即可得到曲顶柱体体积的近似值。柱体体积的近似值。11()()3nniiiibadczznn 表表9-29-2为分别取为分别取10001000、1000010000、5000050000、100000100000、10000001000000得到的近似值，这里要注意的是，由于取点的随
34、得到的近似值，这里要注意的是，由于取点的随机性，对于同样的，每次试验的结果也有可能不相同，机性，对于同样的，每次试验的结果也有可能不相同，但随着的增大，不同的试验得到的近似值的摆动幅度会但随着的增大，不同的试验得到的近似值的摆动幅度会越来越小。越来越小。n n100010001000010000500005000010000010000010000001000000近似值近似值0.31800.31800.31710.31710.31700.31700.31760.31760.31750.3175表表9-29-210:22:56MCM28 现实生活中许多实际问题的求解往往依赖于问题的现实生活中许
35、多实际问题的求解往往依赖于问题的维数，即变量的个数，求解难度随维数的增加呈指数级维数，即变量的个数，求解难度随维数的增加呈指数级增长。有时问题的维数可能高达数百甚至数千，即使使增长。有时问题的维数可能高达数百甚至数千，即使使用当今速度最快的超级计算机，利用传统的数值方法也用当今速度最快的超级计算机，利用传统的数值方法也可能无能为力。蒙特卡罗方法能很好地被用来对付这类可能无能为力。蒙特卡罗方法能很好地被用来对付这类问题，因为此方法的计算复杂性并不依赖于维数。蒙特问题，因为此方法的计算复杂性并不依赖于维数。蒙特卡罗方法还经常被用来检验用其它算法得到的最优解是卡罗方法还经常被用来检验用其它算法得到的
36、最优解是否为真正的最优解，因此蒙特卡罗方法已经成为求解模否为真正的最优解，因此蒙特卡罗方法已经成为求解模型的必备算法之一。型的必备算法之一。10:22:56MCM29 由由n n个不同的字母或数字构成的方阵，每个字母或数个不同的字母或数字构成的方阵，每个字母或数字在每行、每列中都恰好出现一次，这样的方阵被称为字在每行、每列中都恰好出现一次，这样的方阵被称为n n阶拉丁方。如下列两个方阵分别为阶拉丁方。如下列两个方阵分别为3 3阶拉丁方和阶拉丁方和4 4阶拉丁阶拉丁方。方。9.3 9.3 拉丁方与正交设计拉丁方与正交设计123231312123421433412432110:22:56MCM30
37、 在这两个方阵中，每一行均有数字在这两个方阵中，每一行均有数字1 1、2 2、3(3(或或1 1、2 2、3 3、4)4)且均出现一次，每一列同样也均有数字且均出现一次，每一列同样也均有数字1 1、2 2、3(3(或或1 1、2 2、3 3、4)4)且也均出现一次。容易看出，任意且也均出现一次。容易看出，任意n n个不同个不同的字母或数字都至少可以构成一个拉丁方，我们只要在的字母或数字都至少可以构成一个拉丁方，我们只要在第一行按顺序写出第一行按顺序写出n n个字母或数字，然后从第二行起，每个字母或数字，然后从第二行起，每行后退一格，后面的循环转到前面就可以了。行后退一格，后面的循环转到前面就可
38、以了。现我们将两个现我们将两个3 3阶拉丁方阶拉丁方 和和 按照顺按照顺123231312123312231序合并得到一个新的序合并得到一个新的3 3阶方阵阶方阵 ，新方阵中每个，新方阵中每个112233233112321321位置上的元素为原来的两个方阵在对应位置上的元素按位置上的元素为原来的两个方阵在对应位置上的元素按合并次序的得到的数偶。合并次序的得到的数偶。10:22:56MCM31 如果新方阵中的所有元素如果新方阵中的所有元素(数偶数偶)11)11、1212、1313、2121、2222、2323、3131、3232、3333在方阵中都只出现了一次，则我们在方阵中都只出现了一次，则我
39、们称原先的两个称原先的两个3 3阶拉丁方相互正交。阶拉丁方相互正交。一般，若两个一般，若两个n n阶拉丁方按次序合并，两矩阵相应位阶拉丁方按次序合并，两矩阵相应位置的元素置的元素(数字或字母数字或字母)组成的组成的n n2 2个数偶个数偶(或序偶或序偶)在新矩在新矩阵中出现且仅出现一次，则称具有这种性质的拉丁方为阵中出现且仅出现一次，则称具有这种性质的拉丁方为相互正交的拉丁方，简称正交拉丁方，而由两个相互正交的拉丁方，简称正交拉丁方，而由两个n n阶正交阶正交拉丁方按次合并得到的新矩阵称为拉丁方按次合并得到的新矩阵称为n n阶欧拉方阵，简称阶欧拉方阵，简称n n阶欧拉方。阶欧拉方。10:22:
40、56MCM32 并不是所有的拉丁方均具有这种正交性质，如并不是所有的拉丁方均具有这种正交性质，如3 3阶拉阶拉丁方丁方 和和 按序合并之后得到的矩阵按序合并之后得到的矩阵123312231312231123132132321321213213不难验证下列四个不难验证下列四个4 4阶拉丁方中阶拉丁方中 1234214334124321123434124321214312344321214334121234234134124123中，数偶中，数偶1313、2121、3232出现了三次而不是一次，而其余出现了三次而不是一次，而其余6 6个个数偶：数偶：1111，1212，2222，2323，3131
41、，3333均没有出现，因此这两均没有出现，因此这两个拉丁方不正交。个拉丁方不正交。前三个是两两正交的，而最后一个没有和它正交的伙伴。前三个是两两正交的，而最后一个没有和它正交的伙伴。10:22:56MCM33 拉丁方及正交拉丁方的概念在数学及其应用中都很拉丁方及正交拉丁方的概念在数学及其应用中都很有意义。相传有意义。相传1818世纪时，普鲁士国王腓特烈大帝要举行世纪时，普鲁士国王腓特烈大帝要举行一次盛大的阅兵式。当时腓特烈的军队中有一次盛大的阅兵式。当时腓特烈的军队中有6 6个兵种，每个兵种，每个兵种又有个兵种又有6 6种军阶。腓特烈想在前来观礼的各国贵宾面种军阶。腓特烈想在前来观礼的各国贵宾
42、面前炫耀一下，便命令阅兵司令在每一兵种、每一军阶中前炫耀一下，便命令阅兵司令在每一兵种、每一军阶中各挑选一名军官，编成一个由各挑选一名军官，编成一个由3636名军官组成的名军官组成的6 6行行6 6列的列的方阵，要求每一行和每一列中都有各兵种、各军阶的军方阵，要求每一行和每一列中都有各兵种、各军阶的军官。这一要求似乎并不算太高，阅兵司令按要求挑选了官。这一要求似乎并不算太高，阅兵司令按要求挑选了3636名军官，左排右排，可怎么也达不到国王提出的要求。名军官，左排右排，可怎么也达不到国王提出的要求。腓特烈大帝大失所望，一怒之下，革了阅兵司令的职，腓特烈大帝大失所望，一怒之下，革了阅兵司令的职，亲
43、自去指挥排练，但他也始终排不出来。后来，他向大亲自去指挥排练，但他也始终排不出来。后来，他向大数学家欧拉求救，希望他能帮助解决这一困难，这就形数学家欧拉求救，希望他能帮助解决这一困难，这就形成了历史上非常著名的成了历史上非常著名的“3636军官问题军官问题”。10:22:56MCM34 由前面的拉丁方及欧拉方的定义，由前面的拉丁方及欧拉方的定义，“3636军官问题军官问题”实际上是在问：是否存在两个相互正交的实际上是在问：是否存在两个相互正交的6 6阶拉丁方或是阶拉丁方或是否存在否存在6 6阶欧拉方。欧拉一生解决过许多令人望而生畏的阶欧拉方。欧拉一生解决过许多令人望而生畏的数学难题，但这时他竟
44、被这个似乎并不起眼的问题给难数学难题，但这时他竟被这个似乎并不起眼的问题给难住了，很长一段时间没有得到结果。在那个年代，人们住了，很长一段时间没有得到结果。在那个年代，人们已经利用当时的方法找到了已经利用当时的方法找到了3 3阶、阶、4 4阶和阶和5 5阶的正交拉丁方，阶的正交拉丁方，但对于但对于6 6阶正交拉丁方，欧拉苦苦思索，仍然毫无结果。阶正交拉丁方，欧拉苦苦思索，仍然毫无结果。当他发现当他发现2 2阶正交拉丁方也不存在（非常容易验证）时，阶正交拉丁方也不存在（非常容易验证）时，便作了猜想：对任何非负整数便作了猜想：对任何非负整数k k，n=4k+2n=4k+2阶正交拉丁方都阶正交拉丁方
45、都不存在。不存在。10:22:56MCM35 118 118年之后，法国数学家塔瑞严格证明了年之后，法国数学家塔瑞严格证明了6 6阶阶(即即4k+24k+2中中k=1k=1时时)正交拉丁方确实不存在，这也相当于宣告了正交拉丁方确实不存在，这也相当于宣告了“3636军官问题军官问题”的无解。这个结果似乎增加了人们对欧的无解。这个结果似乎增加了人们对欧拉猜想的信任。然而，又过了拉猜想的信任。然而，又过了6262年，数学家帕克却证明年，数学家帕克却证明了了1010阶阶(即即4k+24k+2中中k=2k=2时时)正交拉丁方是存在的，并且将其正交拉丁方是存在的，并且将其构造了出来（如图构造了出来（如图9
46、.89.8所示），这个所示），这个1010阶欧拉方的出现相阶欧拉方的出现相距欧拉提出猜想竟整整过去了距欧拉提出猜想竟整整过去了180180年！这个欧拉方由两个年！这个欧拉方由两个1010阶拉丁方按序合并而成，两个拉丁方均由阶拉丁方按序合并而成，两个拉丁方均由0 0、1 1，9 9这十个数字构成，合并之后得到这十个数字构成，合并之后得到0000，0101，9999共共100100个个数，如果把这些数都看成是整数的话数，如果把这些数都看成是整数的话(00(00看成看成0 0，0909看成看成9 9，其余，其余9090个数偶均看成对应的两位数个数偶均看成对应的两位数)，这，这100100个整个整数的
47、满足上述要求的排列在长达数的满足上述要求的排列在长达180180年的时间内竟被普遍年的时间内竟被普遍认为是不可能的！欧拉的猜想就这样被推翻了。认为是不可能的！欧拉的猜想就这样被推翻了。10:22:56MCM36到了到了19601960年，数学家们彻底解决了欧拉猜想，证明了除年，数学家们彻底解决了欧拉猜想，证明了除了了k=0k=0和和1 1之外，其它的之外，其它的n=4k+2n=4k+2阶正交拉丁方都是存在的。阶正交拉丁方都是存在的。下面就是那个人们经过下面就是那个人们经过180180年努力才找到的由两个年努力才找到的由两个1010阶正阶正交拉丁方构成的矩阵：交拉丁方构成的矩阵：46576870
48、81021324359971943765124029068853932654013819857760421543802709746658923132781689635547910420670579524436908321 1884694133259872105607593022149761084573862811039650873462497500829548762351391764 图图9-89-8 10:22:56MCM37 拉丁方及正交拉丁方的应用非常广泛，下面我们来看拉丁方及正交拉丁方的应用非常广泛，下面我们来看几个简单的实例。几个简单的实例。例例9.39.3 某农场引进了某农场引进了
49、4 4个小麦新品种，要通过较大面积的个小麦新品种，要通过较大面积的栽培试验，看看哪个品种更适合在本地发展。如果把这栽培试验，看看哪个品种更适合在本地发展。如果把这4 4个品种分别种在个品种分别种在4 4块地里，虽然块地里，虽然4 4块地的大气候大致相同，块地的大气候大致相同，但小气候仍然会有差别，很可能会影响到试验结论的准确但小气候仍然会有差别，很可能会影响到试验结论的准确性。此外这性。此外这4 4块地又有不可避免的土壤结构和肥沃程度的块地又有不可避免的土壤结构和肥沃程度的差异，也有可能对结论的准确性产生影响。为此决定将一差异，也有可能对结论的准确性产生影响。为此决定将一块较大面积的试验田划分
50、成块较大面积的试验田划分成4 4行行4 4列共列共1616块，块，4 4种小麦按种小麦按4 4阶阶拉丁方进行种植，使得每种小麦种植在每行每列的一块上，拉丁方进行种植，使得每种小麦种植在每行每列的一块上，且每行每列上均只种植一块。这样处理能做到在进行对比且每行每列上均只种植一块。这样处理能做到在进行对比试验时，气候、土壤结构及肥沃程度等所产生的影响较为试验时，气候、土壤结构及肥沃程度等所产生的影响较为均等，在一定程度上消除了外界环境对实验结果的影响，均等，在一定程度上消除了外界环境对实验结果的影响，使试验结果能更准确地反映出品种的优劣。使试验结果能更准确地反映出品种的优劣。10:22:56MCM
