蒙特卡洛能解决啥

How do you solve an ‘unsolvable’ problem?

您如何解决“无法解决”的问题？

The worlds of data science, mathematical finance, physics, engineering and bioinformatics (amongst many others) readily produce intractable problems. These are problems for which no computationally ‘easy’ solutions are available.

数据科学，数学金融，物理学，工程学和生物信息学(以及许多其他领域)世界容易产生棘手的问题。 这些都是没有计算上的“简便”解决方案的问题。

Luckily, there are methods that can approximate the solutions to these problems with a remarkably simple trick.

幸运的是，有些方法可以用非常简单的技巧来近似解决这些问题。

Monte Carlo methods are a class of methods that can be applied to computationally ‘difficult’ problems to arrive at near-enough accurate answers. The general premise is remarkably simple:

蒙特卡洛方法是一类方法，可以应用于计算上的“困难”问题，以得出几乎足够的准确答案。 总体前提非常简单：

1. Randomly sample input(s) to the problem随机抽样问题的输入
2. For each sample, compute an output对于每个样本，计算输出
3. Aggregate the outputs to approximate the solution汇总输出以近似解决方案

As an analogy, imagine you’re an ant crawling over a large, tiled mosaic. From your vantage point, you have no easy way of working out what the mosaic depicts.

打个比方，假设您是一只蚂蚁，爬在一块大的瓷砖马赛克上。 从有利的角度出发，您没有简便的方法可以计算出马赛克所描绘的内容。

If you started walking around the mosaic and sampling the tiles you visit at random intervals, you’d build up an approximate idea of what the mosaic shows. The more samples you take, the better your approximation will be.

如果您开始在马赛克周围走动并以随机间隔对访问的图块进行采样，那么您将对马赛克显示的内容有一个大概的了解。 您抽取的样本越多，您的近似值越好。

If you could cover every single tile, you’d eventually have a perfect representation of the mosaic. However, this wouldn’t be necessary — after a certain amount of sampling, you’d have a pretty good estimate.

如果您可以覆盖每块瓷砖，那么最终将可以完美呈现马赛克。 但是，这不是必需的-经过一定数量的采样后，您将获得一个不错的估计。

This is exactly how Monte Carlo methods approximate solutions to otherwise ‘unsolvable’ problems.

这正是蒙特卡洛方法近似解决原本“无法解决”的问题的方式。

The name refers to a famous casino in Monaco. It was coined in 1949 by one of the method’s pioneers, Stanislaw Ulam. Ulam’s uncle was reportedly a gambling man, and the connection between the element of ‘chance’ in gambling and in Monte Carlo methods must have been particularly apparent to Stanislaw.

这个名字指的是摩纳哥著名的赌场。 它是由该方法的先驱之一Stanislaw Ulam在1949年发明的。 据报道，乌兰姆的叔叔是一个赌博人，在斯坦尼斯拉夫看来，赌博中的“机会”元素与蒙特卡洛方法之间的联系一定很明显。

The best way to understand a technical concept is to dive right in and see how it works. The rest of this article will show how Monte Carlo methods can solve three interesting problems. The examples will be presented in the Julia programming language.

理解技术概念的最佳方法是直接研究并了解其工作原理。 本文的其余部分将展示蒙特卡洛方法如何解决三个有趣的问题。 这些示例将以Julia编程语言呈现。

Julia介绍 (Introducing Julia)

There are a number of languages you might consider learning if you are interested in specialising in data science. One which has emerged as an increasingly serious option in recent years is a language called Julia.

如果您对数据科学专业感兴趣， 可以考虑学习多种语言 。 近年来，一种日益严重的选择出现了一种叫做Julia的语言。

Julia is a numerical programming language that has seen adoption within a range of quantitative disciplines. It is free to download. There’s also a really neat browser-based interface called JuliaBox, which is powered by Jupyter Notebook.

Julia是一种数字编程语言，已在一系列定量学科中得到采用。 可以免费下载。 还有一个非常整洁的基于浏览器的界面，名为JuliaBox，由Jupyter Notebook支持 。

One of the cool features of Julia we’ll be making use of today is how easily it facilitates parallel computing. This allows you to carry out computations on multiple processes, giving a serious performance boost when done at scale.

我们今天将使用的Julia的很酷的功能之一是它促进并行计算的容易程度。 这样一来，您就可以在多个进程中进行计算，从而可以大规模扩展性能。

平行进行 (Going parallel)

To launch Julia on multiple processes, go to the terminal (or open a new terminal session in JuliaBox) and run the following command:

要在多个进程上启动Julia，请转到终端(或在JuliaBox中打开新的终端会话)并运行以下命令：

$ julia -p 4

This initiates a Julia session on four CPUs. To define a function in Julia, the following syntax is used:

这将在四个CPU上启动Julia会话。 要在Julia中定义函数，请使用以下语法：

function square(x) return x^2 end

That’s right — instead of indentations or curly braces, Julia uses a start-end approach. For-loops are similar:

没错-Julia使用start-end方法来代替缩进或花括号。 For循环类似于：

for i = 1:10 print(i) end

You can of course add whitespace and indentation to aid readability.

您当然可以添加空格和缩进以提高可读性。

Julia’s parallel programming capabilities rely primarily upon two concepts: remote references and remote calls.

Julia的并行编程功能主要依赖于两个概念： 远程引用和远程调用 。

• Remote references are objects that essentially act as named placeholders for objects defined on other processes.远程引用是本质上充当其他进程上定义的对象的命名占位符的对象。
• Remote calls allow processes to call functions on arguments stored on other processes.远程调用允许进程对存储在其他进程上的参数调用函数。

It is important to define functions across all processes. Check out the code below:

在所有流程中定义功能很重要。 查看以下代码：

@everywhere function hello(x)return "Hello " * xendresult = @spawn hello("World!")
print(result)
fetched = fetch(result)
print(fetched)

The @everywhere macro ensures the hello() function is defined across all processes. The @spawn macro is used to wrap a closure around the expression hello("World!"), which is then evaluated remotely on an automatically chosen process.

@everywhere宏可确保在所有进程中定义hello()函数。 @spawn宏用于将闭包包装在表达式hello("World!") ，然后在自动选择的进程上对其进行远程求值。

The result of that expression is instantly returned as a Future remote reference. If you try printing result, you’ll be disappointed. The output of hello("World!") has been evaluated on a different process, and isn’t available here. To make it available, use the fetch() method.

该表达式的结果将立即作为Future远程引用返回。 如果您尝试打印result ，您将感到失望。 hello("World!")的输出已在其他进程上进行了评估，在此处不可用。 要使其可用，请使用fetch()方法。

If spawning and fetching seem like too much to bother with, then you’re in luck. Julia also has a @parallel macro that will take on some of the heavy lifting required for running tasks in parallel.

如果生成和获取似乎太麻烦了，那么您很幸运。 Julia还有一个@parallel宏，它将承担并行运行任务所需的一些繁重工作。

@parallel works either standalone, or with a ‘reducer’ function to collect the results across all processes and reduce them to a final output. Take a look at the code below:

@parallel既可以独立使用，也可以使用“减少器”功能来收集所有过程中的结果并将其减少为最终输出。 看下面的代码：

@parallel (+) for i = 1:1000000000return iend

The for-loop simply returns the value of i at each step. The @parallel macro uses the addition operator as a reducer. It takes each value of i and adds it on to the previous values.

for循环仅在每一步返回i的值。 @parallel宏将加法运算符用作@parallel器。 它获取i每个值，并将其添加到先前的值上。

The result is the sum of the first billion integers.

结果是前十亿个整数的总和。

With that whistle-stop tour of Julia’s parallel programming capabilities in mind, let’s move on to seeing how we can use Monte Carlo methods to solve some interesting example problems.

牢记着Julia的并行编程功能，让我们继续前进，看看如何使用蒙特卡洛方法解决一些有趣的示例问题。

玩彩票 (Playing the lottery)

As a first example, let’s imagine playing a lottery game. The idea is simple — pick six unique numbers between 1 and 50. Each ticket costs, say, £2.

作为第一个示例，让我们想象一下玩彩票游戏。 这个想法很简单-在1至50之间选择六个唯一的数字。每张票的成本为£2。

• If you match all six numbers to those drawn, you win a large prize (£1,000,000)如果您将所有六个数字都与抽奖的数字相匹配，您将赢得一笔大奖(£1,000,000)
• If you match five numbers, you win a medium prize (£100,000)如果您匹配五个数字，您将赢得中奖(100,000英镑)
• If you match four numbers, you win a small prize (£100)如果您匹配四个数字，您将赢得一个小奖(100英镑)
• If you match three numbers, you win a very small prize (£10)如果您匹配三个数字，您将赢得很少的奖金(10英镑)

What would you expect to win if you played this lottery every day for twenty years?

如果您连续二十年每天都玩此彩票，您会期望赢得什么？

You could work this out with pen and paper, by using a little probability theory. But that’d be time consuming! Instead, why not use a Monte Carlo method?

您可以使用一点概率论，用纸和笔解决这个问题。 但这会很耗时！ 相反，为什么不使用蒙特卡洛方法呢？

The approach is almost suspiciously simple — simulate the game over and over many times, and average the outcome.

这种方法几乎是可疑的简单-反复模拟游戏多次，并取平均结果。

Start Julia:

开始Julia：

$ julia -p 4

Now, import the StatsBase package. Use the @everywhere macro to make it available… well, everywhere.

现在，导入StatsBase包。 使用@everywhere宏使其可用……好，随处可见。

using StatsBase@everywhere using StatsBase

Next, define a function that will simulate a single lottery game. The arguments allow you to change the rules of the game, to explore different scenarios.

接下来，定义一个将模拟单个彩票游戏的函数。 这些参数允许您更改游戏规则，以探索不同的场景。

@everywhere function lottery(n, outOf, price)ticket = sample(1:outOf, n, replace = false)draw = sample(1:outOf, n, replace = false)matches = sum(indexin(ticket,draw) .!= 0 )if matches == 6return 1000000 - priceelseif matches == 5return 100000 - priceelseif matches == 4return 100 - priceelseif matches == 3return 10 - priceelsereturn 0 - priceendend

The number of matching numbers is calculated using Julia’s indexin() function. This takes an array, and for each element, returns the index of its position in another array (or zero if the element is not found). Unlike many modern languages, Julia indexes from one, not zero.

匹配数字的数量是使用Julia的indexin()函数计算的。 这需要一个数组，并且对于每个元素，返回其在另一个数组中位置的索引(如果找不到该元素，则返回零)。 与许多现代语言不同，Julia从1开始而不是从0开始索引。

The .!= 0 syntax checks which of these indices are not equal to zero, and returns either true or false for each. Finally, the number of true’s are summed up, giving the total matching numbers.

.!= 0语法检查这些索引中的哪个索引不等于零，并为每个索引返回true或false 。 最后，将true的数量相加，得出匹配的总数。

Now let’s simulate playing the lottery every day for twenty years… ten thousand times in parallel.

现在，让我们模拟二十年来每天玩彩票的过程……并行进行一万次。

winnings = @parallel (+) for i = 1:(365*20*10000)          lottery(6,50,2)endprint(winnings/10000)

Not a great return, huh?

回报不是很大，是吗？

You could extend the code to allow for more advanced rules and scenarios, and see the effect this has on the outcome.

您可以扩展代码以允许使用更高级的规则和方案，并查看其对结果的影响。

Monte Carlo simulations allow for the modelling of considerably more complex situations than this lottery example. However, the approach is much the same as presented here.

与该彩票示例相比，蒙特卡洛模拟可以对更为复杂的情况进行建模。 但是，该方法与此处介绍的方法大致相同。

Let’s see what else Monte Carlo methods allow us to do…

让我们看看蒙特卡洛方法还可以做什么……

pi的值 (The value of pi)

Pi (or π) is a mathematical constant. It is perhaps most famous for its appearance in the formula for the area of a circle:

Pi(或π)是一个数学常数。 它可能以其在圆面积公式中的外观而闻名：

A = πr²

A =πr²

π is an example of an irrational number. Its exact value is impossible to represent as any fraction of two integers. In fact, π is also an example of a transcendental number — there aren’t even any polynomial equations to which it is a solution.

π是无理数的一个例子。 它的确切值不可能表示为两个整数的任何分数。 实际上，π也是一个超越数的示例-甚至没有任何多项式方程可以解决。

You might think this makes obtaining an accurate value for π less-than-straightforward. Or does it?

您可能会认为这使得获得π的准确值小于直接计算。 还是呢？

Actually, you can find a pretty good estimate of π using a Monte Carlo-inspired method. A visual analogy might go as follows:

实际上，您可以使用蒙特卡洛启发的方法找到一个非常好的π估计。 视觉类比可能如下：

• Draw a 2m×2m square on a wall. Inside, draw a circle of radius 1m.在墙上画一个2m×2m的正方形。 在内部，画一个半径为1m的圆。
• Now, take a few steps back and hurl paint randomly at the wall. Count each time the paint lands in the circle.现在，向后走几步，然后在墙上随机喷漆。 计数每次油漆落在圆圈中的时间。
• After a hundred throws, work out what fraction of throws landed in the circle. Multiply this by the area of the square. There is your estimate for π.投掷100次后，算出落在圆圈中的投掷比例。 将此乘以正方形的面积。 有您对π的估计。

The reason this works is very intuitive.When sampling random points from a square containing a circle, the probability of selecting points from within the circle is proportional to the circle’s area.

这样做的原因非常直观：从包含圆的正方形中采样随机点时，从圆内选择点的概率与圆的面​​积成比例。

With enough random samples, we can find a reliable estimate for this proportion, p.

有了足够的随机样本，我们可以找到该比例p的可靠估计。

Now, we know the area of the square is 2×2 = 4m², and we know the area of the circle is π×r². Since the radius r equals 1, the area of the circle is just π.

现在，我们知道了正方形的面积是2×2 =4m²，我们知道了圆的面积是π×R²。 由于半径r等于1，因此圆的面积仅为π。

As we know the area of the square and have an estimate for the proportion p of its area covered by the circle, we can estimate π. Simply multiply p×4.

由于我们知道正方形的面积，并有一个估计值，可以估算出该正方形被圆覆盖的面积的比例p ，因此可以估算出π。 只需乘以p ×4。

The more random samples we throw, the better the estimate p will be. However, the gain in accuracy diminishes as we take more and more samples.

我们抛出的随机样本越多，估计p越好。 但是，随着我们获取越来越多的样本，准确性的增益会降低。

Here’s the Julia code for simulating this example. I ran this in the JuliaBox terminal, using the following command to launch Julia on four CPUs:

这是用于模拟此示例的Julia代码。 我在JuliaBox终端中运行此命令，使用以下命令在四个CPU上启动Julia：

$ julia -p 4

First, define a sampling method.

首先，定义采样方法。

@everywhere function throwPaint(N)hits = 0for i = 1:Nx = rand() ; y = rand()if x^2 + y^2 < 1hits += 1endendreturn float(hits / N * 4)
end

This runs a loop, randomly sampling x and y coordinates between 0 and 1. The if-statement uses the circle equation to check if the points lie within an imaginary circle, counting the number of hits. The function returns the proportion of hits, multiplied by four.

这会运行一个循环，在0和1之间随机采样x和y坐标。if语句使用圆方程式检查点是否位于虚圆内，计算命中次数。 该函数返回命中比例乘以四。

Running this function in parallel will allow for an extremely high number of samples to be drawn, giving much greater precision.

并行运行此功能将允许抽取大量样本，从而提供更高的精度。

Pi = @parallel (+) for i = 1:nworkers()              throwPaint(100000000) / nworkers() endprint(Pi)

The nworkers() method returns the number of CPUs in use (in this case, four). This means each process runs the throwPaint() method, a hundred million times. Overall, this gives us a huge number of samples, and makes for a very precise estimate of the value of π.

nworkers()方法返回正在使用的CPU数(在本例中为四个)。 这意味着每个进程都会运行一亿次throwPaint()方法。 总的来说，这给了我们大量的样本，并且对π值进行了非常精确的估计。

更大的图景：整合 (The bigger picture: Integration)

The estimating π example above is a specific example of a more general use case for Monte Carlo approximation — solving integration problems.

上面的估计π示例是一个更通用的蒙特卡洛近似用例的特定示例-解决积分问题。

Integration is a calculus technique that finds an area defined by a mathematical function. For example, a simple curve might be defined by the function:

积分是一种微积分技术 ，可以找到由数学函数定义的区域。 例如，该函数可能定义了一条简单的曲线：

f(x) → x²

f(x)→x²

And the corresponding graph would be:

相应的图形为：

The area underneath the curve is found by integrating f(x).

通过积分f(x)可以找到曲线下方的面积。

For simpler functions, integration is pretty easy to solve with a little practice. However, for more complicated functions, we need to turn to estimation methods.

对于更简单的功能， 只需一点实践即可很容易地解决集成问题 。 但是，对于更复杂的功能，我们需要转向估算方法。

In low dimensions, the area under a curve can be approximated by relatively simple algorithms, such as the trapezium method.

在低维情况下，可以通过相对简单的算法(例如梯形法)来近似曲线下的面积。

However, because of the curse of dimensionality, this becomes computationally unfeasible in higher dimensions. Monte Carlo-based methods can be used to estimate the area instead.

但是，由于维数的诅咒，在更高的维数中这在计算上变得不可行。 可以使用基于蒙特卡洛的方法来估计面积。

This can be visualized in exactly the same way as the π example above, except the curve need not be defined as a circle. Instead, imagine throwing paint at a unit square containing any arbitrary shape. For example:

可以用与上述π示例完全相同的方式将其可视化，只是不必将曲线定义为圆。 相反，想象一下将油漆投掷到包含任意形状的单位正方形上。 例如：

In higher dimensions, the premise remains the same. The problem is still solved by randomly sampling input values, evaluating them, and aggregating to approximate the solution. Instead of sampling from a circle within a square, imagine sampling a sphere within a cube.

在更高维度上，前提保持不变。 通过随机采样输入值，评估它们并进行汇总以近似解决方案，仍然可以解决该问题。 想象从一个正方形内的一个球体采样，而不是从一个正方形内的一个圆采样。

As a final example, let’s take on a difficult math puzzle.

作为最后一个例子，让我们来解决一个难题。

一个困难的数学难题 (A difficult math puzzle)

Pick two points at random within a unit cube. On average, what is the distance between them?

在一个单位立方体内随机选择两个点。 平均而言，它们之间的距离是多少？

I’ll give you a warning now — the math solution isn’t exactly trivial.

现在，我会警告您- 数学解决方案并不简单 。

However, it is possible to obtain an accurate estimate using — you guessed it — a Monte Carlo method.

但是，可以使用您猜到的蒙特卡洛方法获得准确的估算值。

$ julia -p 4

First, define a sampling method.

首先，定义采样方法。

@everywhere function samplePoints(dimensions)pt1 = []pt2 = []for i = 1:dimensionspt1 = push!(pt1, rand())pt2 = push!(pt2, rand())endreturn [pt1, pt2]end

Now define a function that calculates the distance between the points.

现在定义一个计算点之间距离的函数。

@everywhere function distance(points)pt1 = points[1]pt2 = points[2]arr = []for i = 1:length(pt1)d = (pt2[i] - pt1[i]) ^ 2arr = push!(arr, d)enddist = sqrt(sum(arr))return distend

Finally, run these two functions together in parallel. Instead of reducing to a single output, this time we’ll write each result to a SharedArray object. SharedArray objects allow different processes to access data stored in the same array object.

最后，并行运行这两个功能。 这次我们没有将结果简化为单个结果，而是将每个结果写入SharedArray对象。 SharedArray对象允许不同的进程访问存储在同一数组对象中的数据。

results = SharedArray{Float64}(1000000)
@parallel for i = 1:1000000results[i] = distance(samplePoints(3))endsum(results) / length(results)

You should get an answer very close to 0.6617 — and this is of course the correct answer! By changing the argument passed to samplePoints(), you can solve the generalised problem in however many dimensions you like.

您应该得到非常接近0.6617的答案-这当然是正确的答案！ 通过更改传递给samplePoints()的参数，可以解决您喜欢的多个维度中的广义问题。

接下来是什么？ (What next?)

Hopefully you’ve found this intro to Monte Carlo methods useful!

希望您发现此蒙特卡洛方法介绍很有用！

When implemented correctly, they provide an invaluable tool for data scientists, engineers, financial mathematicians and researchers… and anyone else whose work involves understanding complex systems.

如果正确实施，它们将为数据科学家，工程师，金融数学家和研究人员……以及其他任何需要了解复杂系统的人员提供宝贵的工具。

If you’re interested in learning more about their applications, there’s a ton of resources online. However, the best way to learn is practice! Once you’re comfortable with the basic premise, why not have a go simulating your own Monte Carlo examples?

如果您想了解有关其应用程序的更多信息，则可以在线获取大量资源。 但是，最好的学习方法是练习！ 熟悉基本前提后，为什么不去模拟自己的蒙特卡洛示例呢？

Any feedback or comments, please leave below!

任何反馈或意见，请在下面留下！

翻译自: https://www.freecodecamp.org/news/solve-the-unsolvable-with-monte-carlo-methods-294de03c80cd/

蒙特卡洛能解决啥