如何用python编程实现模拟计算机_Python实现蒙特卡洛模拟（第1部分）【个人翻译自Medium】...-CSDN博客

原文鼓掌数（建议有条件的同学去支持一下原文~）：截至2020/3/4

这是关于学习使用Python进行蒙特卡洛模拟的三部分系列文章的第一个。第一个教程将教您如何进行基本的“粗略”蒙特卡洛，并教您如何使用重要性采样来提高精度。第2部分将介绍臭名昭著的Metropolis算法，第3部分将是专门针对萌芽物理学家的文章（我们将学习如何使用蒙特卡洛模拟来解决量子力学中的问题！）

蒙特卡洛方法是一种广泛使用的启发式技术，可以解决各种常见问题，包括优化和数值积分问题。这些算法通过巧妙地从分布中采样以模拟系统的工作原理来工作。应用范围从解决理论物理问题到预测金融投资趋势。

任务

在本文中，我们将使用Python实现的蒙特卡洛近似法，来以找到以下积分的方案：我们将使用蒙特卡洛方法计算该积分

我借用了这个示例，它是Sunkanta Deb在

您会注意到，该积分无法通过解析求解，因此我们将需要对其进行数值近似。让我们从解决这个问题的简单方法开始：朴素蒙特卡洛。

朴素蒙特卡洛：概念

朴素蒙特卡洛（Crude Monte Carlo）是最容易理解的蒙特卡洛方法，因此我们从这里开始。但是顾名思义，它也是最不准确的。我们稍后会解决此问题。:)

您可能从高中微积分中回忆了以下定理，该定理显示了如何计算函数的平均值（作为证明，请查看此链接）：该定理用于查找函数的平均值（average）

正如我们可以通过积分找到一个函数的平均值一样，我们也可以通过确定一个积分的平均值f来找到该积分的值。蒙特卡洛技术基于此原理：不是估计一个有时候可能无法计算的定积分，而是估计一个定积分的平均值，并以此不断接近定积分的准确值。

而这正是我们要做的！那么我们该怎么做呢？

好吧，最简单的方法是从所有可能的输入值中随机采样输入值x。假设我们有一个简单的线性函数，例如y = 2x，并且我们想要找到y的平均值，其范围为[0,2]。使用10个随机选择的值逼近y = 2x的平均值

为了计算平均值，我们将对所有随机确定的 x求 y ，并将结果取平均值。

这个过程正是朴素蒙特卡洛方法。现在，让我们看看如何在Python中使用此技术解决我们的问题。

朴素蒙特卡洛：应用

你可以找到所有的代码本教程我Github上这里。这是我们代码的开始：

import numpy as np

import math

import random

from matplotlib import pyplot as plt

from IPython.display import clear_output

PI = 3.1415926

e = 2.71828

我们真的不需要太多。numpy将仅用于查找列表的最小元素。math用于定义数学函数，random用于采样，matplotlib将帮助我们直观地观察我们的结果。

一些有用的函数

首先让我们定义一个函数以生成特定范围内的随机数。

def get_rand_number(min_value, max_value):

"""This function gets a random number from a uniform distribution betweenthe two input values [min_value, max_value] inclusivelyArgs:- min_value (float)- max_value (float)Return:- Random number between this range (float)"""

range = max_value - min_value

choice = random.uniform(0,1)

return min_value + range*choice

还有我们的被积函数 f（x）：

def f_of_x(x):

"""This is the main function we want to integrate over.Args:- x (float) : input to function; must be in radiansReturn:- output of function f(x) (float)"""

return (e**(-1*x))/(1+(x-1)**2)

进行朴素蒙特卡洛

实施蒙特卡洛原油计划应该相当简单。我们的算法如下所示：从积分范围中获取随机输入值x

计算被积函数在该点的值

只要您愿意，重复步骤1和2

计算所有这些样本的平均值，然后乘以积分范围

让我们看看对应代码：

def crude_monte_carlo(num_samples=5000):

"""This function performs the Crude Monte Carlo for ourspecific function f(x) on the range x=0 to x=5.Notice that this bound is sufficient because f(x)approaches 0 at around PI.Args:- num_samples (float) : number of samplesReturn:- Crude Monte Carlo estimation (float)"""

lower_bound = 0

upper_bound = 5

sum_of_samples = 0

for i in range(num_samples):

x = get_rand_number(lower_bound, upper_bound)

sum_of_samples += f_of_x(x)

return (upper_bound - lower_bound) * float(sum_of_samples/num_samples)

用N=10000样本进行这种近似计算得出了0.6994。和Wolfram的近似结果的0.696092差的不算太远（我们采取的是朴素法）。

确定我们估计值的方差

但是我们对我们的回答有多自信？我怎么知道10,000个样本足以得到一个很好的近似值。我们可以通过找到估计的方差来量化我们的准确性。方差定义为“与均值的平方距离的平均值”。可以证明它等于以下式子：一般方差方程

方差使我们知道f（x）在x的特定取值范围内的变化幅度。它应该与所使用的样本数无关，但是我们可以通过将方差的平方根除以样本数来计算积分误差。有关方差的更多信息，请查看Adrian Feiguin教授的《蒙特卡洛积分》笔记本，或Wolfram的方差文章。在我们的原始蒙特卡洛模拟中用于计算方差的方程式

这是我们将在仿真中使用的方差方程。让我们看看如何在Python中执行此操作：

def get_crude_MC_variance(num_samples):

"""This function returns the variance fo the Crude Monte Carlo.Note that the inputed number of samples does not neccissarilyneed to correspond to number of samples used in the MonteCarlo Simulation.Args:- num_samples (int)Return:- Variance for Crude Monte Carlo approximation of f(x) (float)"""

int_max = 5 # this is the max of our integration range

# get the average of squares

running_total = 0

for i in range(num_samples):

x = get_rand_number(0, int_max)

running_total += f_of_x(x)**2

sum_of_sqs = running_total*int_max / num_samples

# get square of average

running_total = 0

for i in range(num_samples):

x = get_rand_number(0, int_max)

running_total = f_of_x(x)

sq_ave = (int_max*running_total/num_samples)**2

return sum_of_sqs - sq_ave

蒙特卡洛原油的这种实施给了我一个方差0.266，而误差为0.005。快速思考一下，这效果一点也不差，但是如果我们需要更高的精度怎么办？我们总是可以增加样本数，但是我们的计算时间也会增加。如果，不是使用随机抽样，而是使用更加机灵的抽样方式......这种方法叫做，重要性采样（importance sampling）。

重要性采样：概念

重要性采样是一种在不增加样本数量的情况下减少蒙特卡洛模拟的方差的方法。我们的想法是，与其从整个函数中随机抽样，不如从与函数形状相似的点分布中抽样。

假设您有一个阶跃函数，如下所示：阶跃函数图f（x）及其“权重函数” g（x）

如上所示，我们有一个在[0,2]范围内有效而在[2,6]范围内无效的步进功能。采样10次可能会得出如下估计：f（x）的10个样本估计

这些样本（我发誓是随机的）对应于样本的可能分布，并得出1.8的积分估计值。但是，如果我们取而代之估计我们的函数f(x)和一个特殊的权重函数之间的比值（这个权重函数在任何一个x的取值上都约为f(x) 的取值的一半），结果又会如何呢？如果我们更进一步，将我们的采样更加偏向于那些f(x)的取值波动更加剧烈的x区间（该举措被证明为能够有效减小误差），结果又会如何呢？你将会看到，该比值将会更加接近该积分式的实际取值：2。使用重要性采样估计阶跃函数面积

而用于确定函数g（x）最佳形式的方法则是重要性采样方法。

重要性采样的数学原理

我将简要介绍重要性抽样的数学原理，但是本文的主要目的是编程实现，因此，如果您希望看到更严格的数学过程，请查看我之前提到的Deb文章。

让我们看看是否可以找到这样的g（x）：

基本上，我们希望g（x）看起来像f（x）的一种缩放。

我们还需要g（x）来满足一些条件：g（x）是可积的

g [x ]在[a，b]上为非负数

g（x）的不定积分（我们称为G（x））具有实数逆

[a，b]范围内g（x）的积分必须等于1

在理想情况下，f（x） = k * g（x），其中k为常数。但是，如果f（x） = k * g（x）严格成立，则f（x）将是可积的，若如此我们将无需使用蒙特卡洛模拟方法来估计f（x）的积分，只要使用数学分析方法就能算出该积分。

因此，我们将约定f（x） ≈ k * g（x）。我们当然不会得到一个完美的估计，但是您会发现它的性能比之前的朴素蒙特卡洛估计要好。

我们将如下定义G（x），还将变量更改为r。G(x)的定义变量替换：将自变量从x更改为r（其中r将会在0到1之间的随机抽样）

r将限制在[0,1]范围内。由于g（x）的积分被规定为1，所以G（x）永远不会大于1，因此r永远不会大于1。这很重要，因为稍后在执行模拟时，我们将从[0,1]范围内的r中随机采样。

使用这些定义，我们可以得出以下估计：

这是我们执行蒙特卡洛时要计算的总和。我们将从r中随机取样以产生我们的估计。

简单吧？乍看之下，这没有任何意义，但是请不要相信自己这样过于直觉的想法。我在此有意识地专注于直观，对于数学只是稍微提及。如果您感到困惑，或者想更要严格的数学，请查看我之前提到的资源，直到您相信最终公式。

重要性采样的Python实现

好的，现在我们了解了重要性采样的数学原理，让我们从以前回到问题所在。请记住，我们正在尝试尽可能精确地估计以下积分：我们正在尝试近似的积分式

可视化我们的问题

让我们首先为g（x）权重函数生成一个模板。我想可视化我的函数f（x），所以我们将使用matplotlib来实现：

xs = [float(i/50) for i in range(int(50*PI*2))]

ys = [f_of_x(x) for x in xs]

plt.plot(xs,ys)

plt.title("f(x)");绘制此函数将有助于我们确定g（x）的最佳结构

好了，我们可以看到，我们的函数在区间[0,3-ish]上极其积极，而在区间 [3-ish, inf] 上极其不积极。因此，让我们看看是否可以找到参数化的函数模板，以模仿这种性质。Deb的建议是使用以下函数：建议的权重函数模板[1]

找到A和λ的理想值后，我们将能够得到f（x）和最佳权重函数g（x）的图像：

您会看到g（x）在许多方面不能理想地模仿f（x）的形状。但是这种程度其实是可以接受的。和积分函数粗略相似的g（x）仍然可以有效地降低您估计的方差。您也可以随意尝试其他形式的权重函数g（x），看看是否可以找到更好的解决方案。

确定 g（x）中的参数

在我们可以进行模拟之前，需要找到参数 λ和A 的最佳取值。我们可以使用对g（x）的归一化限制来找到A（λ）：g（x）的归一化条件可以用来证明A（λ） = λ

现在，我们需要做的就是找到理想的λ，然后就能得到理想的g（x）。

为此，让我们以0.5步长在范围[0.05,3.0]上计算不同λ的方差，并使用方差最低的λ 。当使用重要性抽样时，我们计算f（x）与g（x）之比的方差。

然后，我们将要使用此方程式来近似计算积分：

我们将多次计算不同λ的方差，并相应地更改权重函数。之后，我们将使用最优λ计算具有最小方差的积分。

该算法将如下所示：从λ= 0.05开始

计算方差

增加λ

重复步骤2和3，直到达到最后一个λ

选择方差最小的λ，这是您的最佳λ

使用具有此λ的重要性抽样蒙特卡洛计算积分

在代码中，找到最佳λ看起来像这样：

# this is the template of our weight function g(x)

def g_of_x(x, A, lamda):

e = 2.71828

return A*math.pow(e, -1*lamda*x)

def inverse_G_of_r(r, lamda):

return (-1 * math.log(float(r)))/lamda

def get_IS_variance(lamda, num_samples):

"""This function calculates the variance if a Monte Carlousing importance sampling.Args:- lamda (float) : lamdba value of g(x) being testedReturn:- Variance"""

A = lamda

int_max = 5

# get sum of squares

running_total = 0

for i in range(num_samples):

x = get_rand_number(0, int_max)

running_total += (f_of_x(x)/g_of_x(x, A, lamda))**2

sum_of_sqs = running_total / num_samples

# get squared average

running_total = 0

for i in range(num_samples):

x = get_rand_number(0, int_max)

running_total += f_of_x(x)/g_of_x(x, A, lamda)

sq_ave = (running_total/num_samples)**2

return sum_of_sqs - sq_ave

# get variance as a function of lambda by testing many

# different lambdas

test_lamdas = [i*0.05 for i in range(1, 61)]

variances = []

for i, lamda in enumerate(test_lamdas):

print(f"lambda {i+1}/{len(test_lamdas)}: {lamda}")

A = lamda

variances.append(get_IS_variance(lamda, 10000))

clear_output(wait=True)

optimal_lamda = test_lamdas[np.argmin(np.asarray(variances))]

IS_variance = variances[np.argmin(np.asarray(variances))]

print(f"Optimal Lambda: {optimal_lamda}")

print(f"Optimal Variance: {IS_variance}")

print(f"Error: {(IS_variance/10000)**0.5}")

您会看到，使用10,000个样本运行此优化代码会产生λ值1.65，其方差0.0465，误差0.022，哇！使用重要性抽样可以使我们在样本数量相同的情况下将误差减少2倍。

开始模拟

现在，我们要做的就是使用优化的g（x）函数运行仿真，过程一切顺利。这是代码的样子：

def importance_sampling_MC(lamda, num_samples):

A = lamda

running_total = 0

for i in range(num_samples):

r = get_rand_number(0,1)

running_total += f_of_x(inverse_G_of_r(r, lamda=lamda))/g_of_x(inverse_G_of_r(r, lamda=lamda), A, lamda)

approximation = float(running_total/num_samples)

return approximation

# run simulation

num_samples = 10000

approx = importance_sampling_MC(optimal_lamda, num_samples)

variance = get_IS_variance(optimal_lamda, num_samples)

error = (variance/num_samples)**0.5

# display results

print(f"Importance Sampling Approximation: {approx}")

print(f"Variance: {variance}")

print(f"Error: {error}")

运行这段代码使我得到的近似值0.6983与Wolfram提供的伟大答案0.696更加接近。

总结

在本教程中，我们学习了如何执行蒙特卡洛模拟来估计定积分。我们同时使用了粗略方法和重要性抽样方法，发现重要性抽样显着降低了误差。

稍后再查看本教程的第2部分，我们将介绍如何在Monte Carlo模拟中使用Metropolis算法！

祝你好运，并且一如既往地随意提问和评论:)

参考资料和进一步阅读

建议阅读的其他资料：时雨：重要性采样（Importance Sampling）zhuanlan.zhihu.com忘荃：50行代码实现抽象画【蒙特卡洛抽样】zhuanlan.zhihu.com通常蒙特卡罗方法可以粗略地分成两类：一类是所求解的问题本身具有内在的随机性，借助计算机的运算能力可以直接模拟这种随机的过程。例如在核物理研究中，分析中子在反应堆中的传输过程。中子与原子核作用受到量子力学规律的制约，人们只能知道它们相互作用发生的概率，却无法准确获得中子与原子核作用时的位置以及裂变产生的新中子的行进速率和方向。科学家依据其概率进行随机抽样得到裂变位置、速度和方向，这样模拟大量中子的行为后，经过统计就能获得中子传输的范围，作为反应堆设计的依据。