随机算法
随机算法分为两大类:蒙特卡罗算法和拉斯维加斯算法,都是以著名的赌城命名的,且都是通过随机采样尽可能找到最优解。
这两类随机算法之间的选择,往往受到问题的局限。
如果问题要求在有限采样内,必须给出一个解,但不要求是最优解,那就要用蒙特卡罗算法。
反之,如果问题要求必须给出最优解,但对采样没有限制,那就要用拉斯维加斯算法。
比如,机器围棋程序,因为每一步棋的运算时间、堆栈空间都是有限的,而且不要求最优解,所以ZEN涉及的随机算法,肯定是蒙特卡罗式的。
一、蒙特·卡罗算法
(一)背景
此算法诞生的背景是:
1、曼哈顿计划,有极大的计算需求。
2、计算机刚开始发展,最适合做计算。
蒙特卡洛算法理论基础是概率论,实际就是暴力计算逼近理想结果。正是在以上两个背景下,它刚好得到了极大的应用和发展。
(二)概念
蒙特卡洛是一个地名,位于赌城摩纳哥,象征概率。蒙特卡洛(Monte Carlo)方法是由大名鼎鼎的数学家冯·诺伊曼提出的,诞生于上世纪40年代美国的“曼哈顿计划”。原理是通过大量随机样本,去了解一个系统,进而得到所要计算的值。蒙特·卡罗方法在金融工程学,宏观经济学,计算物理学(如粒子输运计算、量子热力学计算、空气动力学计算)等领域应用广泛。
(三)实例1:蒙特卡罗法计算圆周率
1、计算圆周率π的原理
- 一个圆半径R,它有一个外切正方形边长2R。
可以知道:
- 圆面积=Pi R2
- 正方形面积 2R x 2R=4R2
- 从这个正方形内随机抽取一个点,对这个点的要求是在正方形内任意一点的概率平均分布。
那么这个点在圆以内的概率大概就是pi*R2/4R2=pi/4 - 生成若干个这样的点,利用平面上两点间距离公式计算这个点到圆心的距离来判断是否在圆内。设两个点A、B以及坐标分别为
- ,
则A和B两点之间的距离为:
当我们使用足够多的点来进行统计时,我们得到的概率值十分接近pi/4
这样就可以得到pi值
java实现代码:
public class Test {
public static void main(String[] args) {
int n = 100000;//投点次数
int count = 0;//投点成功的次数
for(int i=0;i<n;i++){
double x = Math.random();
double y = Math.random();
//计算点到圆心的距离
if(Math.sqrt((Math.pow(x,2)+Math.pow(y,2)))<1.0){
count++;
}
if(count%3000==0){ //每成功3000次输出一次
System.out.println(4.0*count/n);
}
}
}
}
Python实现代码:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
n = 10000 # 投点次数
r = 1.0 # 圆的半径
a,b = (0.,0.) #圆心
# 正方形区域
x_min, x_max = a - r, a + r
y_min, y_max = b - r, b + r
# 在正方形区域内随机投点
x = np.random.uniform(x_min, x_max, n) #均匀分布
y = np.random.uniform(y_min, y_max, n)
#计算点到圆心的距离
d = np.sqrt((x-a)**2 + (y-b)**2)
#统计落在圆内点的数目
res = sum(np.where(d < r, 1, 0))
#计算pi的近似值(Monte Carlo:用统计值去近似真实值)
pi = 4 * res / n
print('pi: ',pi)
#画个图
fig = plt.figure()
axes = fig.add_subplot(111)
axes.plot(x, y,'ro',markersize = 1)
plt.axis('equal') #防止图像变形
circle = Circle(xy=(a,b), radius=r ,alpha=0.5)
axes.add_patch(circle)
plt.show()
(四)蒙特卡罗法求定积分的近似值
1、在一个1×1的正方形区域里,使用蒙特卡洛方法,随机在这个正方形里面产生大量随机点(数量为n),计算有多少点(数量为res)落在函数下方区域内,res/n就是所要求的积分值,也即红色区域的面积。
2、Python实现代码:
import numpy as np
import matplotlib.pyplot as plt
def f(x):
#return x**2 + 4*x*np.sin(x)
return x**3/3.0+4.0*np.sin(x) - 4.0*x*np.cos(x)
n = 30000
#矩形边界区域
x_min, x_max = 0.0, 1.0
y_min, y_max = 0.0, 1.0
#在矩形区域内随机投点x
x = np.random.uniform(x_min, x_max, n)#均匀分布
y = np.random.uniform(y_min, y_max, n)
#统计落在函数y=x^2下方的点
res = sum(np.where(y < f(x), 1 ,0))
#计算定积分的近似值
integral = res / n
print('integeral: ', integral)
# 画个图
fig = plt.figure()
axes = fig.add_subplot(111)
axes.plot(x, y,'ro',markersize = 1)
plt.axis('equal') # 防止图像变形
axes.plot(np.linspace(x_min, x_max, 10), f(np.linspace(x_min, x_max, 10)), 'b-')
# 显示函数图像
plt.show()