蒙特卡罗算法

蒙特卡洛法又称统计模拟法,得名于摩纳哥著名的赌场蒙特卡洛。该方法解决问题的基本步骤是:

  • 构造一个概率模型,使得要解决的问题正好是概率模型的一个参数,比如概率模型的期望值
  • 根据构建的概率模型生成样本
  • 从样本中构建一个估计器并给出问题的近似解

下面以pi和自然常数的解为例,简单介绍一下模型卡罗法

π%5Cpi

  • 构建概率模型
  • 给定一个边长为2的正方形和期内切圆,内切圆和正常性的面积比率$ \frac{内切圆面积}{正方形面积} = \frac{\pi * 1^2}{2 * 2} = \frac{\pi}{4}$
  • 如果在正方形内部随机产生一定数量的点,那么一个点落入内接圆的概率等于该内接圆与正方形面积的比值
  • X%2CY%20%5Csim%20U%28-1.0%2C%201.0%29,然后P%28X%5E2%20%2B%20Y%5E2%20%5Cle%201%29%20%3D%20%5Cfrac%7B%5Cpi%7D%7B4%7D
  • 采样
  • X%20%5Csim%20U%28-1.0%2C%201.0%29%2C%20Y%20%5Csim%20U%28-1.0%2C%201.0%29,构成样本%28x_1%2Cy_1%29%2C...%2C%28x_n%2Cy_n%29,其中n为样本量
  • 构建估算器
  • 样本落入圆圈的概率P%28X%5E2%20%2B%20Y%5E2%20%3C%201%29%20%3D%20%5Cfrac%7B%5Cpi%7D%7B4%7D%5Capprox%20%5Cfrac%7B%5Csum%201%5Cquad%20if%20%5Cquad%20x_i%5E2%20%2B%20y%5E2%20%5Cle%201%20%7D%7Bn%7D
  • %5Cpi%20%5Capprox%204%20%2A%20%5Cfrac%7B%20%5Csum_i1%5Cquad%20if%5Cquad%20x_i%5E2%20%2B%20y_i%5E2%20%3C%3D%201%20%7D%7Bn%7D%20%3D%204%20%2A%20%5Cfrac%7Bm%7D%7Bn%7D
import matplotlib.pyplot as plt
import numpy as np

n = 10000

x = np.random.uniform(-1.0, 1.0, n)
y = np.random.uniform(-1.0, 1.0, n)


m = 0
inner_x = []
inner_y = []
outer_x = []
outer_y = []
for i in range(n):
    if x[i]*x[i] + y[i]*y[i] <= 1:
        m += 1
        inner_x.append(x[i])
        inner_y.append(y[i])
    else:
        outer_x.append(x[i])
        outer_y.append(y[i])
        
plt.figure(figsize=(8, 8))
plt.scatter(inner_x, inner_y, color='red', s = 10)
plt.scatter(outer_x, outer_y, color = 'blue', s = 10)


pi = 4 * m / n
print('n=%d, pi = %f' % (n, pi))


请添加图片描述
请添加图片描述

自然常数 e

请添加图片描述

  • 构建概率模型
  • %5Cint_1%5E2%20%5Cfrac%7B1%7D%7Bx%7Ddx%20%3D%20ln2%20-%20ln1%20%3D%20ln2
  • %5Cint_1%5E2%5Cfrac%7B1%7D%7Bx%7Ddx是上图中阴影部分的面积,它的面积是用蒙特卡洛法计算的
  • 要在正方形中产生均匀分布的随机点,点落入阴影部分的概率等于阴影面积与矩形面积的比p%20%3D%20%5Cfrac%7BA%E7%9A%84%E9%9D%A2%E7%A7%AF%7D%7B%E7%9F%A9%E5%BD%A2%E9%9D%A2%E7%A7%AF%7D%20%3D%20%5Cfrac%7B%5Cint_1%5E2%5Cfrac%7B1%7D%7Bx%7Ddx%7D%7B%282-1%29%2A%201%7D
  • 采样
  • X%20%5Csim%20U%281.0%2C%202.0%29%2C%20Y%20%5Csim%20I%280.0%2C%201.0%29,抽取样本%28x_1%2Cy_1%29%2C...%2C%28x_n%2Cy_n%29, n为样本容量
  • 构建估算器
  • %5Chat%7Bp%7D%20%3D%20%5Cfrac%7B%5Csum_i%201%5Cquad%20if%5Cquad%20y_i%20%5Cle%20%5Cfrac%7B1%7D%7Bx_i%7D%7D%7Bn%7D%20%3D%20%5Cint_1%5E2%5Cfrac%7B1%7D%7Bx%7Ddx%20%3D%20ln2%20%3D%20%5Cfrac%7Bm%7D%7Bn%7D
  • ln2%20%3D%20%5Cfrac%7Blog_22%7D%7Blog_2e%7D%20%3D%20%5Cfrac%7Bm%7D%7Bn%7D
  • e%3D2%5E%7B%5Cfrac%7Bn%7D%7Bm%7D%7D
import matplotlib.pyplot as plt
import numpy as np


n = 1000
x = np.random.uniform(1.0, 2.0, n)
y = np.random.uniform(0.0, 1.0, n)

m = 0
inner_x = []
inner_y = []
outer_x = []
outer_y = []
for i in range(n):
    if y[i] <= 1 / x[i]:
        m += 1
        inner_x.append(x[i])
        inner_y.append(y[i])
    else:
        outer_x.append(x[i])
        outer_y.append(y[i])
        
plt.figure(figsize=(10, 10))
plt.scatter(inner_x, inner_y, color='red', s = 10)
plt.scatter(outer_x, outer_y, color = 'blue', s = 10)

print("n=%d, e = %f" % (n, np.power(2, n / m)))


请添加图片描述
请添加图片描述

版权声明:本文为博主jony0917原创文章,版权归属原作者,如果侵权,请联系我们删除!

原文链接:https://blog.csdn.net/gaofeipaopaotang/article/details/123031918

共计人评分,平均

到目前为止还没有投票!成为第一位评论此文章。

(0)
乘风的头像乘风管理团队
上一篇 2022年2月23日 上午11:29
下一篇 2022年2月23日 上午11:50

相关推荐