蒙特卡罗(Monte Carlo)方法,或称计算机随机模拟方法,是一种基于“随机数”的计算方法。这一方法源于美国在第一次世界大战进研制原子弹的“曼哈顿计划”。该计划的主持人之一、数学家冯·诺伊曼用驰名世界的赌城—摩纳哥的Monte Carlo—来命名这种方法,为它蒙上了一层神秘色彩。
' O. R1 l5 @; R' q+ ^9 Y* e% E Monte Carlo方法的基本思想很早以前就被人们所发现和利用。早在17世纪,人们就知道用事件发生的“频率”来决定事件的“概率”。19世纪人们用投针试验的方法来决定圆周率π。本世纪40年代电子计算机的出现,特别是近年来高速电子计算机的出现,使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能。
0 c! M6 S% v( j6 {9 j* D* ^ 考虑平面上的一个边长为1的正方形及其内部的一个形状不规则的“图形”,如何求出这个“图形”的面积呢?Monte Carlo方法是这样一种“随机化”的方法:向该正方形“随机地”投掷N个点落于“图形”内,则该“图形”的面积近似为M/N。 2 ]3 ? x8 ^, b- {
可用民意测验来作一个不严格的比喻。民意测验的人不是征询每一个登记选民的意见,而是通过对选民进行小规模的抽样调查来确定可能的优胜者。其基本思想是一样的。
5 u7 D6 D0 F" R: {, g% v* C v# o& I, T 科技计算中的问题比这要复杂得多。比如金融衍生产品(期权、期货、掉期等)的定价及交易风险估算,问题的维数(即变量的个数)可能高达数百甚至数千。对这类问题,难度随维数的增加呈指数增长,这就是所谓的“维数的灾难”(Course Dimensionality),传统的数值方法难以对付(即使使用速度最快的计算机)。Monte Carlo方法能很好地用来对付维数的灾难,因为该方法的计算复杂性不再依赖于维数。以前那些本来是无法计算的问题现在也能够计算量。为提高方法的效率,科学家们提出了许多所谓的“方差缩减”技巧。
! u! |' I6 W8 k" x; V" t 另一类形式与Monte Carlo方法相似,但理论基础不同的方法—“拟蒙特卡罗方法”
. x) i3 ?2 M5 e(Quasi-Monte Carlo方法)—近年来也获得迅速发展。我国数学家华罗庚、王元提出的“华—王”方法即是其中的一例。这种方法的基本思想是“用确定性的超均匀分布序列(数学上称为Low Discrepancy Sequences)代替Monte Carlo方法中的随机数序列。对某些问题该方法的实际速度一般可比Monte Carlo方法提出高数百倍,并可计算精确度。 具体实现的matlab代码: b8 T( c/ _1 ? Q6 D5 M2 }
---------------------------------------------------------------------------------------------------6 B9 Y. c: }: t& O2 f% c0 t
function val = ballvol(n, m)5 O# \, d' e! ^% z; r) F! u, {
% BALLVOL Compute volume of unit ball in R^n* s: z) V# u) H
%
5 S; h; z n$ z% Computes the volume of the n-dimensional unit ball 7 W0 S/ ~2 }% z+ |& }
% using monte-carlo method.
; ?# }2 G. g ~0 v) P% usage: val = BallVol(n, m)
0 l5 o: C) C+ ~9 p" P9 ~% where: n = dimension * m& l' I( K: _1 H6 b- e
% m = number of realisations
1 S3 E/ J3 }4 F# y9 Z" L* P% If the second argument is omitted, 1e4 is taken as default for m.
9 o& Q' p1 R: b- D+ S2 y- C; l4 j' f- _. s$ S( U
% (c) 1998, Rolf Krause, krause@math.fu-berlin.de \ a- p& R5 a% F, |+ Q
, h! G/ S. ?. ~# t) a9 r3 B- f" R
M = 1e4;: W: E* T! {8 ]4 `8 C
error = 0;+ d0 f9 e7 }7 U
if(nargin <1 | nargin > 2), error('wrong number of arguments'); end
9 ^2 G4 t3 c8 Nif nargin == 2, M = m; end ( W! s5 E1 q; M
. t) ^& {6 d5 C( BR = rand(n, M);$ N7 R+ D2 d' C4 U! U" `
in = 0;
/ H% `% k" [' ]for i=1:M
% C# J( p1 M( W+ Rif(norm(R(:,i),2) <= 1.0), in = in+1; end
$ A0 g% g5 G+ r1 ?: ], I& send2 u& ]1 R- w9 g9 E
! A, M( D P1 }( L: K& }- P9 wval = 2^n*in/M;0 B( ]1 E6 T$ b+ u& ~9 H
--------------------------------------------------------------------------------------------------
9 X1 O- R. h: U# H8 s+ Y |