蒙特卡洛算法

Code Writers

于 2023-09-02 17:03:21 发布

阅读量2.8k

点赞数 4

文章标签：算法 python java 开发语言概率论 eclipse

本文链接：https://blog.csdn.net/m0_71905144/article/details/132640296

版权

文章目录

一蒙特卡罗算法简介
二求圆周率 $π$ 的python实例演示
三算法的优缺点
四算法的应用
总结

一蒙特卡罗算法简介

蒙特·卡罗方法（Monte Carlo method），也称统计模拟方法，它是一种思想或者方法的统称，而不是严格意义上的算法。蒙特卡罗方法的起源是1777年由法国数学家布丰（Comte de Buffon）提出的用投针实验方法求圆周率（具体算法见文末的好文推荐），在20世纪40年代中期，由于计算机的发明结合概率统计理论的指导，从而正式总结为一种数值计算方法，其主要是用随机数来估算计算问题。

蒙特卡罗算法基本思想

投针实验这个例子有点晦涩，换个简单的图文并茂的例子–求圆周率。计算圆周率时可以考虑将一个单位圆放在一个正方形中，从而将求解圆周率转化为计算出圆和正方形面积的比例。

蒙特卡罗方法的基本思想是这样：假想你有一袋豆子，把豆子均匀地朝这个正方形上撒，然后数落在圆内的豆子数占正方形内豆子数的比例，即可计算出圆的面积，近而计算出 $π$ 。而且当豆子越小，撒的越多的时候，结果就越精确。蒙特卡罗的理论依据是概率论中的大数定律

蒙特卡罗算法的基本步骤

蒙特卡罗算法一般分为三个步骤，包括构造随机的概率的过程，从构造随机概率分布中抽样，求解估计量。

1 构造随机的概率过程

对于本身就具有随机性质的问题，要正确描述和模拟这个概率过程。对于本来不是随机性质的确定性问题，比如计算定积分，就必须事先构造一个人为的概率过程了。它的某些参数正好是所要求问题的解，即要将不具有随机性质的问题转化为随机性质的问题。如本例中求圆周率的问题，是一个确定性的问题，需要事先构造一个概率过程，将其转化为随机性问题，即豆子落在圆内的概率，而π就是所要求的解。

2 从已知概率分布抽样

由于各种概率模型都可以看作是由各种各样的概率分布构成的，因此产生已知概率分布的随机变量，就成为实现蒙特卡罗方法模拟实验的基本手段。如本例中采用的就是最简单、最基本的（0，1）上的均匀分布，而随机数是我们实现蒙特卡罗模拟的基本工具。

3 求解估计量

实现模拟实验后，要确定一个随机变量，作为所要求问题的解，即无偏估计。建立估计量，相当于对实验结果进行考察，从而得到问题的解。如求出的近似π就认为是一种无偏估计。

估计读者读完这三个带有严重概率统计学的段落依旧一头雾水，不着急，结合python实例就明白了。

二求圆周率 $π$ 的python实例演示

借助计算机程序可以生成大量均匀分布的坐标点，接着统计出图形内的点数，通过它们占总点数的比例和坐标点生成范围的面积就可以求出的近似值。

import random
import math
import getopt
import sys
import matplotlib.pyplot as plt

def main():
    total = 0
    print('Start experiment: ')
    n = 0

 #Repeat try the estimate pi, until break it down manually
 while True:
        n=int(input(sys.argv))
 for i in range(n):
            x = random.uniform(-1,1)
            y = random.uniform(-1,1)
 #x^2 + y^2 <=1, means (x, y) in the cycle
 if math.sqrt(x ** 2 + y ** 2) <= 1.0:
                total += 1
                plt.plot(x, y, 'ro')
 # x^2 + y^2 > 1, means (x, y) out of the cycle
 else:
                plt.plot(x, y, 'b*')
        mypi = 4.0 * total / n
 #plot the x, y and close it after 5 seconds
        plt.ion()
        plt.pause(5)
        plt.close()

        print('Iteration Times = ', n, 'PI estimate value = ', mypi)
        print('math.pi = ', math.pi)
        print('Errors = ', abs(math.pi - mypi) / math.pi)