0 建议学时

4学时

1 引入

1.1 随机数与采样

客观世界的某些行为，结果具有随机性：
掷骰子、投硬币；
等待公交车的时间；
种子发芽的比例；
…

1.2 随机数函数

1.2.1 random模块

Python的random模块中提供了若干生成随机数的函数
如random.random()可以产生区间[0,1)内的一个随机数

import random
random.random()  0.2865887089067386
random.random()  0.613263367620591
random.random()  0.04800896336544469

注意：计算机只能产生’‘伪随机数’'，这些随机数本身仍然是由一个确定的算法生成的。

关于random库

random.random() 产生均匀分布于[0, 1)的随机数
random.uniform(a,b) 产生均匀分布于[a, b)的随机数

# 生成100个[-1, 1)的随机数
import random
N = 100
x = range(N)
y = [random.uniform(-1, 1) for item in x]

for语句简写

[x for x in range(10)]
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9][x+3 for x in range(10)]
[3, 4, 5, 6, 7, 8, 9, 10, 11, 12]

随机函数的使用

import matplotlib.pyplot as plt
import random
N = 100
x = range(N)
y = [random.uniform(-1, 1) for item in x]
plt.plot(x, y, '+')
plt.show()

在这里插入图片描述

1.2.2 批量产生随机数

random库中的random函数每次只能产生单个随机数
使用numpy库的random模块可以一次性生成大量的随机数

from numpy import randomrandom.random()  #0.37452288986512927
array = random.random(size=5)
print(array) #array([0.93312885, 0.70160214, 0.9727459, 0.1363383,0.54146683])random.uniform(-1,1)  #0.8232171988008607
array = random.uniform(-1,1,size=5)
print(array) #array([ 0.94529425,-0.6128988,0.29927762,0.65928358,0.11272987])

1.3 Numpy

Numpy：Numerical Python，Numpy是高性能科学计算和数据分析的基础包，其部分功能如下：
– ndarray，一个具有矢量算术运算和复杂广播能力的快速且节省空间的多维数组。
– 用于对整组数据进行快速运算的标准数学函数（无需编写循环）。
– 用于读写磁盘数据的工具以及用于操作内存映射文件的工具。
– 线性代数、随机数生成以及傅里叶变换功能。
– 用于集成由C、C++、Fortran等语言编写的代码的工具。

在这里插入图片描述

1.3.1 安装Numpy

在这里插入图片描述

1.3.2 Numpy库

在这里插入图片描述

1.3.3 arange(): 创建数组

import numpy as nparray = np.arange(5) #左闭右开[0,5)
print(array) #[0 1 2 3 4]
array = np.arange(1, 5) #左闭右开[1,5)
print(array) #[1 2 3 4]
array = np.arange(0, 1, 0.1)
print(array) #[0.  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9]

1.3.4 zeros和ones

zeros和ones分别可以创建指定长度或者形状的全0或全1数组

array = np.zeros((3, 2))
[[0. 0.][0. 0.][0. 0.]]array = np.ones((3, 2))
print(array)
[[1. 1.][1. 1.][1. 1.]]array = np.eye(2)
print(array)
[[1. 0.][0. 1.]]

1.3.5 Numpy运算

import numpy as np 
x=np.array([1,2,3,4,5]) #[1,2,3,4,5]
x=x*2                   #[2,4,6,8,10]
print(x>5)              #[False,False,True,True,True]

在这里插入图片描述

1.3.6 Numpy库速度对比

import numpy as np
import timemy_arr = np.arange(1000000)
my_list = list(range(1000000))
beginTime = time.perf_counter()for _ in range(10):my_arr2 = my_arr * 2tNp = time.perf_counter() -beginTimebeginTime = time.perf_counter()
for _ in range(10):my_list2 = [x * 2 for x in my_list]
tList = time.perf_counter() -beginTimeprint("tNp = ", tNp)
print("tList = ", tList)# tNp =  0.009103899999999998 (9ms)
# tList =  0.3303001 (330ms)

1.3.7 批量产生随机数

[a,b]之间的随机整数
random和numpy都提供了整型随机数函数

import random
r=random.randint(1,100) #[1,100]import numpy as np
r=np.random.randint(1,100,10) #[1,100)# array([33, 53, 24, 12, 49, 24, 63, 74, 29, 94])
r=np.random.random_integers(1,100,10) #[1,100]# array([21, 32, 54, 17, 47, 84, 52, 76, 100, 73])

示例：统计掷骰子出现某点数的频率
掷一枚骰子会随机出现1-6点，掷10000次，出现6的概率是多少？

import random
N = 10000
lt = [random.randint(1, 6) for i in range(N)]
M = lt.count(6)
print('6点%d次，共掷了%d次' % (M, N))
print('出现6点的概率:', M/N)

使用numpy

import numpy as np  
N = 10000
lt = np.random.randint(1, 7, N)
six = [i for i in lt if i == 6]
print('6点%d次，共掷了%d次' % (len(six), N))
print('出现6点的概率:', len(six)/N)

掷骰子（使用随机数组）：使用数组比较和numpy中的数组求和更高效

import numpy as np
N = 10000
eyes = np.random.randint(1, 7, N)
success = eyes == 6
M = np.sum(success)
print('6点%d次，共掷了%d次' % (M, N))
print('出现6点的概率:', M/N)

1.3.8 随机数程序的调试

随机数程序的调试困难，因为每次执行会生成不同的数据序列
可以设置随机种子，使随机数的生成序列固定random.seed(10)
缺省情况下，系统以当前时间作为种子

import random
random.seed(10)
['%.3f' % random.random() for i in range(6)]
#['0.571', '0.429', '0.578', '0.206', '0.813', '0.824']random.seed(10)
['%.3f' % random.random() for i in range(6)]
#['0.571', '0.429', '0.578', '0.206', '0.813', '0.824']

2 蒙特卡洛(Monte Carlo)方法

2.1 问题引入

如何求不规则图形的面积？
在这里插入图片描述
应用：蒙特卡洛剂量计算被广泛视为放疗领域，剂量计算金标准。通过随机模拟大量粒子与物质相互作用的全过程，蒙特卡洛方法准确计算，反映了真实的剂量分布，可以实现5分钟内快速准确的蒙特卡洛剂量计算，剂量计算的统计偏差小于2%。
在这里插入图片描述
常见应用：

估算某事件发生的概率；
估算不规则区域的面积/体积；
蒙特卡洛积分。

2.2 方法介绍

基本思想：通过随机抽样的方法，以随机事件出现的频率估计其概率，或者以抽样的数字特征估算随机变量的数字特征，并将其作为问题的近似解。
在这里插入图片描述
$¼的圆面积=(圆内点数/总点数)∗正方形面积\color{red}{¼的圆面积 = (圆内点数/总点数)*正方形面积}$

2.2 应用1：计算不规则区域的面积

函数 $𝑒^{x^𝟐}$ 在区间 $[a, b] (a > 0)$ 上的曲线如下图所示，计算曲线下与 $x$ 轴围成的阴影部分面积
在这里插入图片描述
【思路】

在矩形区域 $B=[1,2]×[𝒆^{a^𝟐},𝒆^{b^𝟐}]$ 内随机采样 $N$ 个点；
计算落在阴影区域的点 (满足 $0<=y<= 𝒆^{𝒙^𝟐}$ ) 的数目(设为 $M$ 个)；
阴影部分的面积为 $M / N \times a re a (B)$ 。

【练习1】求阴影部分的面积：
$(x−1)2+(y−1)2≤1x2+y2>22\begin{array}{l} (x-1)^{2}+(y-1)^{2} \leq 1 \\ x^{2}+y^{2}>2^{2} \end{array}$
在这里插入图片描述
【练习2】计算不规则区域的体积
求坐标在原点的球
$x2+y2+z2≤4x^2 + y^2 + z^2 \leq 4$
与圆柱
$(x−1)2+y2≤1(x-1)^2 + y^2\leq 1$
的相交部分（称为Viviani体）的体积。
【分析】
在哪个区域采样？
如何判断采样点落在相交区域？
【思路1】相交区域内的点/球外切正方体内的点

import numpy
B = 2; N = 1000000; M = 0X = numpy.random.uniform(-B, B, size=N)
Y = numpy.random.uniform(-B, B, size=N)
Z = numpy.random.uniform(-B, B, size=N)c1=(X**2+Y**2+Z**2<B**2)
c2=((X-1)**2+Y**2<=(B/2)**2)
#c3=numpy.int16(c1)+numpy.int16(c2)
#M=numpy.count_nonzero(c3==2)c3=c1 & c2
M=numpy.sum(c3)V = (2*B)**3
print(M/N*V)

【思路2】相交区域内的点/球内点

from numpy import random
B = 2; N = 10000;
X = random.uniform(-B, B, size=N)
Y = random.uniform(-B, B, size=N)
Z = random.uniform(-B, B, size=N)
N_inS=0; N_inC=0
for i in range(N):if X[i]**2 + Y[i]**2 + Z[i]**2 <= B**2:N_inS += 1if (X[i]-1)**2 + Y[i]**2 <= (B/2)**2:N_inC += 1
V = 4*pi/3*(2**3)
print(N_inC/N_inS*V)

2.3 应用2：蒙特卡洛积分

2.3.1 问题背景

对于区间[a,b]上的可积函数f，可以使用牛顿-莱布尼兹公式：
$∫abf(x)dx=F(x)∣ab\int_{a}^{b} f(x) \mathrm{d} x=\left.F(x)\right|_{a} ^{b}$
其中 $F$ 是其原函数。
然而，某些情形下 $F (x)$ 并不存在解析表达式，如当 $f(x)=x^x$ , $f(x)=e^{x^2}$ ,或者 $f(x)=sin⁡(x)xf(x)=\frac{\sin(x)}{x}$ 等。
可尝试使用蒙特卡洛方法近似计算定积分的值

2.3.2 蒙特卡洛积分

定积分：G是曲线 $y = f (x)$ 与 $x$ 轴之间在 $\in [a,b]$ 上的几何图形, 则G的面积值就是 $∫abf(x)dx\int_{a}^{b} f(x) \mathrm{d} x$ , 设 $B = [a, b] \times [0, m]$ , 则有
在这里插入图片描述
$∫abf(x)dx≈MNm(b−a)\int_{a}^{b} f(x) \mathrm{d} x \approx \frac{M}{N} m(b-a)$
其中， $N$ 为采样点数， $M$ 为落在曲线与x轴之间区域的点的数目。
【例子1】

import numpy as np
def MCint_area(f, a, b, n, fmax):below = 0x = np.random.uniform(a, b,n)y = np.random.uniform(0, fmax,n)yCurve=f(x)below=np.sum(y<yCurve)area = below/n*(b-a)*fmaxreturn areadef f1(x):return np.e**(x**2)area=MCint_area(lambda x:np.e**(x**2), 1, 2, 100000, np.e**(2**2))
print(area)area=MCint_area(f1, 1, 2, 100000, np.e**(2**2))
print(area)

附录 Numpy运算

逻辑运算

np.logical_and([True, False], [False, False])      &
np.logical_or([True, False], [False, False])       |
np.logical_not([True, False])

比较运算

np.greater([4,2],[2,2])
np.greater_equal([4, 2], [2, 2])      
np.less([1, 2], [2, 2])
np.less_equal([4, 2, 1], [2, 2, 2])
np.equal([1.,2], [1., 3])
np.not_equal([1.,2], [1., 3])