蒙特卡罗算法：模拟

最常用的模拟技术是蒙特卡罗分析，这种分析为每个计划活动确定一种活动持续时间概率分布，然后利用这些分布计算出整个项目持续时间可能结果的概率分布。—— 来自PMP考点

结合昨晚看到的一篇文章《蒙特卡罗算法到底是什么》进行学习。

蒙特卡罗算法到底是什么？

常见于金融行业，或者其他领域。

它可以模拟出很多场景，并且模拟出来的数据，相对精准。

pi值的模拟

pi 是一个无理数，无限不循环，早在南北朝时期，我国数学家祖之冲得出精确到小数点后7位的结果。

模拟代码：

import random as random

first_count = 100000
count_s = 0
j = 0

while j < first_count:
    x = random.uniform(-1,1)
    y = random.uniform(0,1)
    if x**2+y**2 < 1:
        count_s = count_s + 1
    j = j + 1

print(count_s/j*4)

## 3.14228

在模拟超过10w次之后，基本已经趋于稳定，3.14。

具体解释：向1个正方形和内切圆投点，看落在圆内的概率为多少？

积分

计算区域内面积，利用蒙特卡罗求函数的积分。

模拟代码：

import random

count = 0
count_s = 0
x1 = 0

for i in range(0, 10000):
    x = random.uniform(-1,1)
    y = random.uniform(0, 1)
    count = count + 1
    if y < -(x**2) + 1:
        count_s = count_s + 1
    elif y == -(x**2) + 1:
        x1 = x1 + 1
print((count_s+x1/2)/count*2)

## 1.3368

总结

通过蒙特卡罗模拟，生成一系列符合预期要求的随机数，就可以模拟出一个十分接近实际值的近似值，一般适用于对数值计算精度要求不高的场景。

学习代码

# -*- coding: utf-8 -*-
"""
Created on Mon Nov 18 17:31:07 2019

@author: Hider
"""

import numpy as np
import math
import random
from time import time
import matplotlib.pyplot as plt

t0 = time()
x1 = 0
x2 = 1
count = 0
count_s = 0
list_x = []
list_y = []

## 计算积分
for i in range(1, 101):
    list_x.append(i * 10)
print(list_x[0])

for j in (list_x):
    for i in range(0, j):
        x = random.uniform(-1, 1)
        y = random.uniform(0, 1)
        count = count + 1
        if y < -(x**2) + 1:
            count_s = count_s + 1
        elif y == -(x**2) + 1:
            x1 = x1 + 1
    list_y.append((count_s+x1/2)/count*2)
print("time = ", time() - t0)
print("S = ", list_y[10])


## 绘图
print(len(list_y)) # 100
print(len(list_x)) # 100

plt.rcParams['savefig.dpi'] = 300 # 图片像素
plt.rcParams['figure.dpi'] = 300 # 分辨率
plt.plot(list_x, list_y)
plt.hlines(y = 1.33333, xmin = -1, xmax = max(list_x), colors = "c", linestyles = "dashed")
plt.show()

## 计算pi值
list_p = []
list_p_x = []
max_count = 100000000
first_count = 10
rate = 2
count_s = 0
j = 0

while first_count < max_count:
    print(first_count)
    while j < first_count:
        x = random.uniform(-1, 1)
        y = random.uniform(0, 1)
        if x**2 + y**2 < 1:
            count_s = count_s + 1
        j = j + 1
    list_p_x.append(first_count)
    list_p.append(count_s/first_count*4)
    j = 0
    count_s = 0
    first_count = first_count * rate
print("count_s = ", count_s)
print("pi = ", list_p)
print("pi x = ", list_p_x)

plt.xlim(0, first_count/2)
plt.ylim(3, 3.3)
plt.plot(list_p_x, list_p)
plt.hlines(y = np.pi, xmin = first_count, xmax = list_p_x, colors = "c", linestyles = "dashed")
plt.show()