蒙特卡罗模拟 - python实现

系统 1362 0

基本原理

蒙特卡罗方法通过抓住事物运动的几何数量和几何特征,利用数字方法来加以模拟,即进行一种数字模拟实验。它是以一个概率模型为基础,按照这个模型所描绘的过程,通过模拟实验的结果,作为问题的近似解。

主要步骤如下:
1.构造或描述概率过程
2.实现从已知概率分布抽样
3.建立各种估计量

示例一:π值的计算

            
              
                import
              
               numpy 
              
                as
              
               np

              
                import
              
               pandas 
              
                as
              
               pd

              
                import
              
               matplotlib
              
                .
              
              pyplot 
              
                as
              
               plt

              
                %
              
               matplotlib inline


              
                # π的计算
              
              

n
              
                =
              
              
                30000
              
              

r
              
                =
              
              
                1.0
              
              
                #半径
              
              
a
              
                ,
              
              b
              
                =
              
              
                (
              
              
                0.0
              
              
                ,
              
              
                0.0
              
              
                )
              
              
                #圆心
              
              
xmin
              
                ,
              
              xmax
              
                =
              
              a
              
                -
              
              r
              
                ,
              
              a
              
                +
              
              r
ymin
              
                ,
              
              ymax
              
                =
              
              b
              
                -
              
              r
              
                ,
              
              b
              
                +
              
              r

x
              
                =
              
              np
              
                .
              
              random
              
                .
              
              uniform
              
                (
              
              xmin
              
                ,
              
              xmax
              
                ,
              
              n
              
                )
              
              
y
              
                =
              
              np
              
                .
              
              random
              
                .
              
              uniform
              
                (
              
              ymin
              
                ,
              
              ymax
              
                ,
              
              n
              
                )
              
              

fig
              
                =
              
              plt
              
                .
              
              figure
              
                (
              
              figsize
              
                =
              
              
                (
              
              
                6
              
              
                ,
              
              
                6
              
              
                )
              
              
                )
              
              
axes
              
                =
              
              fig
              
                .
              
              add_subplot
              
                (
              
              
                1
              
              
                ,
              
              
                1
              
              
                ,
              
              
                1
              
              
                )
              
              
plt
              
                .
              
              plot
              
                (
              
              x
              
                ,
              
              y
              
                ,
              
              
                'ro'
              
              
                ,
              
              markersize
              
                =
              
              
                1
              
              
                )
              
              
plt
              
                .
              
              axis
              
                (
              
              
                'equal'
              
              
                )
              
              

d
              
                =
              
              np
              
                .
              
              sqrt
              
                (
              
              
                (
              
              x
              
                -
              
              a
              
                )
              
              
                **
              
              
                2
              
              
                +
              
              
                (
              
              y
              
                -
              
              b
              
                )
              
              
                **
              
              
                2
              
              
                )
              
              
res
              
                =
              
              
                sum
              
              
                (
              
              np
              
                .
              
              where
              
                (
              
              d
              
                <
              
              r
              
                ,
              
              
                1
              
              
                ,
              
              
                0
              
              
                )
              
              
                )
              
              
                print
              
              
                (
              
              
                '落在圆内的点有%i个'
              
              
                %
              
               res
              
                )
              
              

pi 
              
                =
              
              
                4
              
              
                *
              
              res
              
                /
              
              n

              
                print
              
              
                (
              
              
                "π的近似值为:"
              
              
                ,
              
              pi
              
                )
              
              
                #绘制圆形
              
              
                from
              
               matplotlib
              
                .
              
              patches 
              
                import
              
               Circle
circle
              
                =
              
              Circle
              
                (
              
              xy
              
                =
              
              
                (
              
              a
              
                ,
              
              b
              
                )
              
              
                ,
              
              radius
              
                =
              
              r
              
                ,
              
              alpha
              
                =
              
              
                0.5
              
              
                ,
              
              color
              
                =
              
              
                'r'
              
              
                )
              
              
axes
              
                .
              
              add_patch
              
                (
              
              circle
              
                )
              
              
plt
              
                .
              
              grid
              
                (
              
              
                True
              
              
                ,
              
              linestyle
              
                =
              
              
                '--'
              
              
                ,
              
              linewidth
              
                =
              
              
                '0.5'
              
              
                )
              
              
plt
              
                .
              
              show
              
                (
              
              
                )
              
            
          

示例二:计算积分y=x^2

            
              n
              
                =
              
              
                10000
              
              
                #投点次数
              
              
                #矩形区域边界
              
              
x_min
              
                ,
              
              x_max
              
                =
              
              
                0.0
              
              
                ,
              
              
                1.0
              
              
y_min
              
                ,
              
              y_max
              
                =
              
              
                0.0
              
              
                ,
              
              
                1.0
              
              
                #在矩形区域内随机投点
              
              
x
              
                =
              
              np
              
                .
              
              random
              
                .
              
              uniform
              
                (
              
              x_min
              
                ,
              
              x_max
              
                ,
              
              n
              
                )
              
              
y
              
                =
              
              np
              
                .
              
              random
              
                .
              
              uniform
              
                (
              
              y_min
              
                ,
              
              y_max
              
                ,
              
              n
              
                )
              
              
                #创建函数y=x**2
              
              
                def
              
              
                f
              
              
                (
              
              x
              
                )
              
              
                :
              
              
                return
              
               x
              
                **
              
              
                2
              
              
                #统计落在函数y=x^2图像下方的点的数目
              
              
res
              
                =
              
              
                sum
              
              
                (
              
              np
              
                .
              
              where
              
                (
              
              y
              
                <
              
              f
              
                (
              
              x
              
                )
              
              
                ,
              
              
                1
              
              
                ,
              
              
                0
              
              
                )
              
              
                )
              
              
                #计算定积分的近似值
              
              
integral
              
                =
              
              res
              
                /
              
              n

              
                print
              
              
                (
              
              
                'integral:'
              
              
                ,
              
              integral
              
                )
              
              
                #绘制散点图
              
              
fig
              
                =
              
              plt
              
                .
              
              figure
              
                (
              
              figsize
              
                =
              
              
                (
              
              
                6
              
              
                ,
              
              
                6
              
              
                )
              
              
                )
              
              
axes
              
                =
              
              fig
              
                .
              
              add_subplot
              
                (
              
              
                1
              
              
                ,
              
              
                1
              
              
                ,
              
              
                1
              
              
                )
              
              
axes
              
                .
              
              plot
              
                (
              
              x
              
                ,
              
              y
              
                ,
              
              
                'ro'
              
              
                ,
              
              markersize
              
                =
              
              
                1
              
              
                )
              
              
plt
              
                .
              
              axis
              
                (
              
              
                'equal'
              
              
                )
              
              
                #绘制y=x^2面积图
              
              
xi
              
                =
              
              np
              
                .
              
              linspace
              
                (
              
              
                0
              
              
                ,
              
              
                1
              
              
                ,
              
              
                100
              
              
                )
              
              
yi
              
                =
              
              xi
              
                **
              
              
                2
              
              
plt
              
                .
              
              plot
              
                (
              
              xi
              
                ,
              
              yi
              
                ,
              
              
                '--k'
              
              
                )
              
              
plt
              
                .
              
              fill_between
              
                (
              
              xi
              
                ,
              
              yi
              
                ,
              
              
                0
              
              
                ,
              
              color
              
                =
              
              
                'r'
              
              
                ,
              
              alpha
              
                =
              
              
                0.5
              
              
                ,
              
              label
              
                =
              
              
                'area'
              
              
                )
              
              
plt
              
                .
              
              grid
              
                (
              
              
                )
              
            
          

蒙特卡罗模拟 - python实现_第1张图片


更多文章、技术交流、商务合作、联系博主

微信扫码或搜索:z360901061

微信扫一扫加我为好友

QQ号联系: 360901061

您的支持是博主写作最大的动力,如果您喜欢我的文章,感觉我的文章对您有帮助,请用微信扫描下面二维码支持博主2元、5元、10元、20元等您想捐的金额吧,狠狠点击下面给点支持吧,站长非常感激您!手机微信长按不能支付解决办法:请将微信支付二维码保存到相册,切换到微信,然后点击微信右上角扫一扫功能,选择支付二维码完成支付。

【本文对您有帮助就好】

您的支持是博主写作最大的动力,如果您喜欢我的文章,感觉我的文章对您有帮助,请用微信扫描上面二维码支持博主2元、5元、10元、自定义金额等您想捐的金额吧,站长会非常 感谢您的哦!!!

发表我的评论
最新评论 总共0条评论