首页 > 编程语言 >随机模拟——蒙特卡洛算法的Python实现

随机模拟——蒙特卡洛算法的Python实现

时间:2023-12-13 13:22:47浏览次数:92  
标签:Python random 问题 算法 随机 doors 蒙特卡洛 模拟

蒙特卡洛方法是一类基于随机抽样的数值计算技术,通过模拟随机事件的概率过程,从而近似计算复杂问题的数学期望或积分。其核心思想是通过大量的随机抽样来逼近问题的解,从而在随机性中获得问题的统计特性。蒙特卡洛方法广泛应用于概率统计、物理学、金融工程、生物学等领域。
在蒙特卡洛模拟中,通过生成符合特定分布的随机样本,计算这些样本的函数值,然后利用样本均值或频率等统计量来估计问题的解。这种方法的优势在于对于高维、复杂的问题,蒙特卡洛方法通常能够提供更为准确的数值结果。随机模拟是蒙特卡洛方法的一种具体应用,它通过对问题中的不确定性因素进行随机抽样,从而模拟问题的多种可能性。在金融风险评估、气象预测、粒子物理实验等领域,随机模拟的蒙特卡洛方法为复杂问题的求解提供了一种强大而灵活的工具。

一、随机模拟的导入

模拟是指把某一现实的或抽象的系统的某种特征或部分状态,用另一系统(称为模拟模型)来代替或近似。为了解决某问题,把它变成一个概率模型的求解问题,然后产生符合模型的大量随机数,对产生的随机数进行分析从而求解问题,这种方法叫做随机模拟方法,又称为蒙特卡洛(MonteCarlo)方法。
例如,一个交通路口需要找到一种最优的控制红绿灯信号的办法,使得通过路口的汽车耽搁的平均时间最短,而行人等候过路的时间不超过某一给定的心理极限值。十字路口的信号共有四个方向,每个方向又分直行、左转、右转。因为汽车和行人的到来是随机的,我们要用随机过程来描述四个方向的汽车到来和路口的行人到来过程。理论建模分析很难解决这个最优化问题。但是,我们可以采集汽车和行人到来的频率,用随机模拟方法模拟汽车和行人到来的过程,并模拟各种控制方案,记录不同方案造成的等待时间,通过统计比较找出最优的控制方案。
随机模拟中的随机性可能来自模型本身的随机变量,比如上面描述的汽车和行人到来,也可能是把非随机的问题转换为概率模型的特征量估计问题从而用随机模拟方法解决。(用随机模拟估计圆周率) 为了计算圆周率\(\pi\)的近似值可以使用随机模拟方法。
如果向正方形\(D = \{(x,y): x \in [-1,1], y \in [-1,1]\}\)内随机等可能投点,落入单位圆\(C = \{(x,y): x^2 + y^2 \leq 1 \}\)的概率为面积之比\(p =\frac{\pi}{4}\)。如果独立重复地投了\(N\)个点,落入\(C\)中的点的个数\(\xi\)的平均值为\(E\xi = pN\),由概率的频率解释,

\[\begin{aligned} \frac{\xi}{N} \approx& \frac{\pi}{4}, \quad \pi \approx \hat \pi = \frac{4 \xi}{N} \end{aligned} \]

可以这样给出\pi的近似值。

# 蒙特卡罗法(统计试验法)
import random # 导入随机模块
S = 1e6 # 变量S为试验总次数(值设置得越大,PI的计算越准确,即频率越逼近于概率)
N = 0 # 变量N用于统计落在圆内的试验点的个数
for i in range(int(S)):
    x = random.random() # 获取0-1之间的随机数
    y = random.random() # 获取0-1之间的随机数
    d = (x-0.5)**2+(y-0.5)**2 # 计算试验点到圆心的欧式距离的平方
    if d<=0.5**2: # 通过比较试验点到圆心的欧式距离与圆半径的大小,判断该点是否在圆内
        N+=1
    else:
        pass
PI = 4*N/S
print(PI)

二、Python的随机数生成函数

random模块是Python标准库中的一个模块,用于生成各种类型的随机数。它包含了许多函数和方法,可以用于生成伪随机数、从序列中获取随机元素、洗牌等功能。

随机函数 说明
random.random() 返回随机生成的一个浮点数,范围在[0,1)之间
random.uniform(a, b) 返回随机生成的一个浮点数,范围在[a, b)之间
np.random.rand(d0, d1, …, dn) 返回一个或一组浮点数,范围在[0, 1)之间
np.random.normal(loc=a, scale=b, size=()) 返回满足条件为均值=a, 标准差=b的正态分布(高斯分布)的概率密度随机数
np.random.randn(d0, d1, … dn) 返回标准正态分布(均值=0,标准差=1)的概率密度随机数
np.random.standard_normal(size=()) 返回标准正态分布(均值=0,标准差=1)的概率密度随机数
random.sample(k) 从总体序列或集合中随机选取k个唯一的元素
np.random.randint(a, b, size=(), dtype=int) 返回在范围在[a, b)中的随机整数(含有重复值)
random.randrange(a, b, step=c) 在指定范围内,在指定的基数和步长值形成的集合中获取1个随机数值
random.choice(x) 从指定的序列x中随机获取一个数据
random.shuffle(x) 将一个列表x中的元素打乱,随机排序(俗称:洗牌)
random.seed() 设定随机种子

三、蒙特卡洛算法的基本工作步骤

蒙特卡罗算法一般分为三个步骤,分别为构造随机的概率的过程,从构造随机概率分布中抽样,求解估计量。
构造随机的概率过程:对于本身就具有随机性质的问题,要正确描述和模拟这个概率过程。对于本来不是随机性质的确定性问题,比如计算定积分,就必须事先构造一个人为的概率过程了。它的某些参数正好是所要求问题的解,即要将不具有随机性质的问题转化为随机性质的问题。如本例中求圆周率的问题,是一个确定性的问题,需要事先构造一个概率过程,将其转化为随机性问题,即豆子落在圆内的概率,而π就是所要求的解。
从已知概率分布抽样:由于各种概率模型都可以看作是由各种各样的概率分布构成的,因此产生已知概率分布的随机变量,就成为实现蒙特卡罗方法模拟实验的基本手段。如本例中采用的就是最简单、最基本的(0,1)上的均匀分布,而随机数是我们实现蒙特卡罗模拟的基本工具。
求解估计量:实现模拟实验后,要确定一个随机变量,作为所要求问题的解,即无偏估计。建立估计量,相当于对实验结果进行考察,从而得到问题的解。如求出的近似\(\pi\)就认为是一种无偏估计。

四、随机模拟的应用

4.1计算定积分

利用python计算函数\(y =  x^2\)在[0,1]区间的定积分。

import random
 
m = 1000000
n = 0
for i in range(m):
    x = random.random()
    y = random.random()
    if y >= x ** 2:
        n += 1
r = n / m
print("r = {}".format(r))

4.2三门问题

三门问题(Monty Hall problem)亦称为蒙提霍尔问题、蒙特霍问题或蒙提霍尔悖论,大致出自美国的电视游戏节目Let's Make a Deal。问题名字来自该节目的主持人蒙提·霍尔(Monty Hall)。参赛者会看见三扇关闭了的门,其中一扇的后面有一辆汽车,选中后面有车的那扇门可赢得该汽车,另外两扇门后面则各藏有一只山羊。当参赛者选定了一扇门,但未去开启它的时候,节目主持人开启剩下两扇门的其中一扇,露出其中一只山羊。主持人其后会问参赛者要不要换另一扇仍然关上的门。问题是:换另一扇门会否增加参赛者赢得汽车的机率?如果严格按照上述的条件,即主持人清楚地知道,自己打开的那扇门后是羊,那么答案是会。不换门的话,赢得汽车的几率是1/3。换门的话,赢得汽车的几率是2/3。虽然该问题的答案在逻辑上并不自相矛盾,但十分违反直觉。

def three_doors_question(times: int):
    import random
    not_change = 0
    change = 0
    for i in range(0, times):
        doors = [1, 2, 3]  # 1是汽车
        random.shuffle(doors)
        first_choice = doors[random.randint(0, 2)]
        doors.remove(first_choice)
        # 如果汽车在参与者选后剩下的里面,由主持人排除一个错误答案,剩下汽车
        if 1 in doors:
            doors = [1]
        # 如果汽车已经被参与者选了,主持人随机排除剩下一个错误答案
        else:
            doors.remove(random.choice((2, 3)))
        if first_choice == 1:
            not_change = not_change + 1
        if doors[0] == 1:
            change = change + 1
 
    print(
        f'模拟次数:{times}, 不更改选择赢得汽车的概率为:{not_change / times}, '
        f'更改选择赢得汽车的概率为:{change / times}')
 
 
three_doors_question(10000)
three_doors_question(100000)
three_doors_question(500000)
three_doors_question(1000000)

4.3 随机过程的模拟

几何布朗运动的随机微分方程如下,这意味着我们是在等价鞅测度下进行操作:

\[dS_t=\mu S_tdt+\sigma S_tdW_t \]

这里\(W_t\)是一个维纳过程,或者说是布朗运动,而漂移百分比\(\mu\)和波动百分比\(\sigma\)则为常量。
对于上述随机微分方程进行欧拉离散化可以得到离散时间模型:

\[S_t=S_{t-\Delta t } \exp^{(\mu-{\frac{1}{2}} \sigma^2)\Delta t+\sigma \epsilon_t \sqrt{\Delta t}} \]

其中\(S_t\):证券在\(t\)时刻的价格;\(S_{t-\Delta t}\):证券在\(t-\Delta t\)时刻的价格;\(\mu\):证券收益率的期望值;\(\sigma\):证券收益率的波动率;\(\epsilon_t\):服从正态分布(期望为0,方差为1)。
模拟证券初始价格为200(日收益率均值为0.003,波动率为0.25),时间为1年,步长以日为单位,次数为10000次的几何布朗运动价格。

import matplotlib.pyplot as plt
from pylab import mpl
mpl.rcParams['font.sans-serif']=['SimHei']
mpl.rcParams['axes.unicode_minus']=False
import numpy as np
import numpy.random as npr

S0=100  #初始价格
u=0.03  #收益率的均值
sigma=0.24  #收益率的波动率
T=1  #年化时间长度
I=10000  #模拟的次数
M=252   #年化时间分段数目
dt=T/M  #模拟的每小步步长
S=np.zeros((M+1,I))
S[0]=S0
for t in range(1,M+1):
    S[t]=S[t-1]*np.exp((u-0.5*sigma**2)*dt+sigma*np.sqrt(dt)*npr.standard_normal(I))

#绘制直方图
plt.hist(S[-1],bins=50)
plt.title("模拟最终价格的频数图",fontsize=13)
plt.xlabel(u'最终价格',fontsize=13)
plt.ylabel(u'频数',fontsize=13)
plt.xticks(fontsize=13)
plt.yticks(fontsize=13)
plt.grid(True)
plt.show()

总结

随机模拟,尤其是蒙特卡洛方法,能够通过模拟多种潜在的市场情景,提供更全面、鲁棒的需求预测。以金融衍生品定价为例,蒙特卡洛方法能够在面对复杂金融工具时提供更为精确的估值。在传统的金融市场中,期权定价等问题往往难以用解析方法求解,因为这些金融工具的收益受到多个变量的影响,而这些变量往往难以建立简单的数学模型。蒙特卡洛方法通过模拟这些变量的随机变动,从而估计金融工具的价格。具体而言,假设我们要估计某个股票期权的价格,该期权的价值受到股票价格、波动率、无风险利率等多个因素的影响。我们可以使用蒙特卡洛方法模拟这些因素的未来走势,随机生成多个未来可能的市场情景,计算每种情景下期权的收益,然后取平均值作为期权的估值。这样一来,我们就能够考虑到各种不确定性因素,更全面地评估金融工具的风险。

参考文献

  1. 【Python】蒙特卡洛模拟
  2. 10 随机模拟介绍
  3. Python实现常见随机过程的模拟

标签:Python,random,问题,算法,随机,doors,蒙特卡洛,模拟
From: https://www.cnblogs.com/haohai9309/p/17898787.html

相关文章

  • python——小游戏(ball,bird)
      ball #-*-coding:utf-8-*-"""CreatedonWedDec1309:19:382023@author:kabuqinuo"""importsys#导入sys模块importpygame#导入pygame模块pygame.init()#初始化pygamesize=width,height=640,480#设置窗......
  • 代码随想录算法训练营Day1 | 704.二分查找、27.移除元素
    LeetCode704.二分查找二分查找是一种基础的算法,其核心思想在高中数学中就已经被大家所熟知了,然而对于代码的实现,其细节问题常常令人头疼,比如while循环的条件是什么?middle是该+1还是-1?这些问题需要有一个清晰的认知。题目链接如下:704.二分查找Carl的讲解链接:二分查找算法......
  • 如何在本地跑IPA算法
    参考文章:https://blog.csdn.net/qq_45529538/article/details/1313110971)下载源码https://github.com/ipa320/ipa_coverage_planning 2)安装依赖sudoaptinstallros-noetic-libdlibros-noetic-opengmros-noetic-cob-navigationcoinor-* 3)编译catkin_make-DCATKIN_......
  • Python——第五章:hashlib模块
    hashlib模块hashlib模块是Python中用于加密散列(hash)算法的模块。它提供了对常见的哈希算法(如MD5、SHA-1、SHA-256等)的支持,使得开发者可以轻松地在其应用中进行数据的安全散列。以下是hashlib模块中一些常用的哈希算法:MD5(MessageDigestAlgorithm5):产生128位的哈......
  • 实验 2 处理机调度算法
    1.实验任务1)回顾课本第三章中介绍过的作业或进程调度算法,包括先来先服务、最短作业优先、时间片轮转、多级队列调度和多级反馈队列调度等,介绍上述调度算法的设计原理并分析各自的特点;2)采用高级编程设计语言实现任意一种处理机调度算法;3)下面提供了实现先来先服务调度算法的参考......
  • 【python】文件锁模块fcntl
      #!/usr/bin/python#coding:utf8importosimportsysimporttimeimportfcntl#导入模块classFLOCK(ojbect):def__init__(self,name):""":paramname:文件名"""self.fobj=open(name,'......
  • 【教3妹学编程-算法题】交换得到字典序最小的数组
    3妹:2哥2哥,你有没有看到新闻:周海媚姐因病医治无效,于2023年12月11日离开了我们。2哥 :看到了,真是个悲伤的消息,早晨还看到辟谣,以为没事了呢。3妹:是啊,#再见周芷若#2哥:童年的女神,周海媚演的这版“周芷若”真的很深入人心!被评为“最美周芷若”3妹:哎,人有生老病死,R.I.P.2哥:唉,说点高兴的......
  • Python报错:performance hint: av/logging.pyx:232:5: the GIL to be acquired
     参考:https://stackoverflow.com/questions/77410272/problems-installing-python-av-in-windows-11https://github.com/PyAV-Org/PyAV/issues/1177  ================================  报错信息:C:\Windows.old\Users\chris>pipinstallavDefaultingtouserinstallatio......
  • Python报错:pkg-config could not find libraries ['avformat', 'avcodec', 'av
    参考:https://github.com/PyAV-Org/PyAV/issues/238https://pyav.org/docs/6.1.2/installation.html#mac-os-x  =====================  报错信息:C:\Users\liuxue>pipinstallavCollectingavUsingcachedav-0.3.3.tar.gzInstallingcollectedpackages:avRunning......
  • Python学习多线程、多进程、多协程记录
    一、多线程应用于请求和IO#1.Python中关于使用多线程多进程的库/模块#2.选择并发编程方式(多线程Thread、多进程Process、多协程Coroutine)前置知识: 一、三种有各自的应用场景 1.一个进程中可以启动多个线程 2.一个线程中可以启动多个协程 二、各自优缺点 1......