用蒙特卡罗方法解决“无法解决”问题
你如何解决'无法解决'的问题？
数据科学，数学金融，物理学，工程学和生物信息学等许多领域都容易产生棘手的问题。
幸运的是，有一些方法可以用一个非常简单的技巧来逼近这些问题的解决方案。
蒙特卡洛方法是一类方法，可以应用于计算上的“困难”问题，以达到近乎准确的答案。总的前提非常简单：
随机抽样输入问题
对于每个样本，计算一个输出
汇总输出以近似解决方案
作为一个类比，想象一下你是一只蚂蚁爬在一个大的，拼接的马赛克上。从你的观点来看，你没有简单的方法来描述马赛克所描绘的东西。
如果你开始在马赛克周围走动，并以随机的间隔取样瓷砖，你会建立一个马赛克的大致概念。你采集的样本越多，你的近似值就越好。
如果你能涵盖每一个瓦片，你最终会有一个完美的马赛克的代表。然而，这是不必要的 ，经过 一定量的抽样后，你会有一个很好的估计。
这就是蒙特卡罗方法的近似解，否则无法解决的“问题”。
介绍Julia
Julia是一种数字编程语言，已经在一系列量化学科中得到采用。它可以免费下载。还有一个非常整洁的基于浏览器的界面叫JuliaBox，它由Jupyter Notebook提供支持。
我们今天将要使用的Julia的一个很酷的功能就是简化并行计算。这使您可以在多个进程上执行计算，并在规模完成时提供显着的性能提升。
并行
要在多个进程上启动Julia，请转至终端（或在JuliaBox中打开新的终端会话）并运行以下命令：
$ julia -p 4
这启动了四个CPU上的Julia会话。要在Julia中定义函数，使用以下语法：
function square(x) return x^2 end
Julia使用了一种start-end方法。For-loop类似：
for i = 1:10 print(i) end
您当然可以添加空格和缩进来帮助可读性。
Julia的并行编程能力主要依赖于两个概念：远程引用和远程调用。
远程引用是实质上充当其他进程上定义的对象的命名占位符的对象。
远程调用允许进程调用存储在其他进程上的参数的函数。
在所有流程中定义功能非常重要。查看下面的代码：
@everywhere function hello(x)
return "Hello " * x
end
result = @spawn hello("World!")
print(result)
fetched = fetch(result)
print(fetched)
@everywhere 宏确保在所有进程中定义 hello () 函数。@spawn 宏用于在表达式 hello ( "World! ") 周围环绕一个闭合, 然后在自动选择的进程中远程计算。
该表达式的结果将立即作为Future的远程引用返回。如果你尝试打印result, 你会失望的。您好 ( "World! ") 的输出已在不同的进程中进行了评估, 此处不可用。若要使其可用, 请使用 fetch () 方法。
如果spawning和fetching似乎太困扰了，那么你很幸运。Julia还有一个@parallel宏，它将承担并行运行任务所需的一些繁重工作。
@parallel 可以独立工作, 也可以使用 "reducer" 功能在所有进程中收集结果, 并将其减少到最终输出。查看下面的代码:
@parallel (+) for i = 1:1000000000
return i
end
for循环只是返回i每一步的值。该@parallel宏使用加法运算符作为还原器。它取每个值i并将其添加到以前的值。
玩彩票
作为第一个例子，让我们设想玩彩票游戏。这个想法很简单 - 选择1到50之间的六个唯一数字。每票花费，例如2英镑。
如果您将所有六个数字与所绘制的数字相匹配，您将赢得一笔大奖（£1,000,000）
如果您匹配五个号码，您将赢得中等奖（100,000英镑）
如果你匹配四个号码，你会赢得一个小奖（£100）
如果你匹配三个数字，你会赢得一个非常小的奖金（10英镑）
如果你20年每天玩彩票，你会期望赢得什么？
你可以用笔和纸，通过使用一点概率论来解决这个问题。但那会很耗时！相反，为什么不使用蒙特卡洛方法？
开始julia：
$ julia -p 4
现在，导入StatsBase包。使用@everywhere宏使它可用。
using StatsBase
@everywhere using StatsBase
接下来，定义一个函数来模拟一个彩票游戏。参数允许你改变游戏规则，探索不同的场景。
@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 == 6
return 1000000 - price
elseif matches == 5
return 100000 - price
elseif matches == 4
return 100 - price
elseif matches == 3
return 10 - price
else
return 0 - price
end
end
匹配数字的数量是使用Julia的indexin()函数计算的。这需要一个数组，并且对于每个元素，返回其在另一个数组中的位置的索引（如果未找到该元素，则返回零）。与许多现代语言不同，Julia从1索引，而不是0索引。
该 .!= 0语法检查这些索引中的哪些不等于零，并且返回true或false每个都返回。最后，true总结出数字，给出匹配的总数。
现在让我们模拟20年来每天玩彩票的情况，并行地玩上万次。
winnings = @parallel (+) for i = 1:(365*20*10000)
lottery(6,50,2)
end
print(winnings/10000)
您可以扩展代码以允许使用更高级的规则和场景，并查看其对结果的影响。
蒙特卡洛模拟允许建模比这个彩票例子更复杂的情况。但是，这种方法与此处介绍的方法大致相同。
pi的值
Pi（或π）是一个数学常数。这可能是最出名的，因为它出现在一个圆面积的公式中：
A =πr²
π是一个无理数的例子。它的确切值不可能表示为两个整数的任何部分。事实上，π也是超越数的 一个例子- 甚至没有任何多项式方程可以解决这个问题。
实际上，使用蒙特卡罗启发法可以很好地估计π。视觉比喻可能如下：
在墙上画一个2m×2m的正方形。在里面画一个半径为1m的圆。
现在，退后几步，在墙上随意乱涂油漆。每当油漆落在圆圈内时计数。
经过一百次投球后，计算出投球中投入的比例是多少。乘以正方形的面积。有你对π的估计。
圆的面积大约是正方形的78％
这个工作的原因非常直观。当从包含圆的正方形中抽取随机点时，从圆内选择点的概率与圆的面积成正比。
有了足够的随机样本，我们可以找到这个比例的可靠估计，p。
现在，我们知道了正方形的面积是2×2 =4m²，我们知道了圆的面积是π× [R ²。由于半径r等于1，圆的面积仅为π。
由于我们知道正方形的面积并估计了它所覆盖的区域的比例p，我们可以估计π。简单地乘以p ×4。
我们抛出的随机样本越多，估计p越好。然而，随着我们采取越来越多的样本，准确度的增加会减少。
这里是模拟这个例子的Julia代码。我在JuliaBox终端中运行这个命令，使用以下命令在四个CPU上启动Julia：
$ julia -p 4
首先，定义一种抽样方法。
@everywhere function throwPaint（N）
hits = 0
for i = 1：N
x = rand（）; y = rand（）
if x ^ 2 + y ^ 2 <1
hits + = 1
end
end
return float（hits / N * 4）
end
这将运行一个循环，随机抽样x并y在0和1之间进行坐标。if语句使用圆方程来检查点是否位于一个虚圆内，计算点击次数。该函数返回命中的比例乘以4。
并行运行此功能将允许绘制极高数量的样本，从而提供更高的精度。
Pi = @parallel (+) for i = 1:nworkers()
throwPaint(100000000) / nworkers()
end
print(Pi)
该nworkers()方法返回正在使用的CPU数量（在本例中为四个）。这意味着每个进程运行该throwPaint()方法一亿次。总体而言，这给了我们大量的样本，并且对π的值进行了非常精确的估计。
整合
上面的估计π示例是一个更一般的用于蒙特卡罗逼近的用例的具体示例 - 解决集成问题。
积分是一种微积分技术，可以找到由数学函数定义的区域。例如，一个简单的曲线可能由函数定义：
f（x）→x 2
并且相应的图表是：
通过积分f（x）找到曲线下方的区域。
通过积分f（x）找到曲线下方的区域。
对于更简单的功能来说，整合很容易通过一些练习来解决。但是，对于更复杂的功能，我们需要转向估算方法。
在低维情况下，曲线下面积可以用相对简单的算法来近似，如梯形法。
然而，这在高维方面变得在计算上不可行。蒙特卡罗方法可以用来估计面积。
这可以用与上面的π示例完全相同的方式进行可视化，只是曲线不需要定义为圆。相反，想象在包含任意形状的单位正方形上投掷涂料。例如：
在更高维度上，前提保持不变。问题仍然通过对输入值进行随机抽样，对它们进行评估以及聚合以近似解决方案来解决。设想在立方体内采样一个球体，而不是从一个正方形内的一个圆圈采样。
作为最后一个例子，让我们来看一个困难的数学难题。
一个困难的数学难题
在单位立方体内随机挑选两个点。平均而言，它们之间的距离是多少？
如果你喜欢解决多个积分，那么......对你有好处！对于我们其他人，总是有蒙特卡罗
$ julia -p 4
首先，定义一种抽样方法。
@everywhere function samplePoints（dimensions）
pt1 = []
pt2 = []
for i = 1：dimensions
pt1 = push！（pt1，rand（））
pt2 = push！（pt2，rand（））
end
return [pt1，pt2]
end
现在定义一个函数来计算点之间的距离。
@everywhere function distance(points)
pt1 = points[1]
pt2 = points[2]
arr = []
for i = 1:length(pt1)
d = (pt2[i] - pt1[i]) ^ 2
arr = push!(arr, d)
end
dist = sqrt(sum(arr))
return dist
end
最后，并行运行这两个函数。而不是减少到单个输出，这次我们将每个结果写入一个SharedArray对象。SharedArray对象允许不同的进程访问存储在同一个数组对象中的数据。
results = SharedArray {Float64}（1000000）
@parallel for i = 1：1000000
results [i] = distance（samplePoints（3））
end
sum(results) / length(results)
你应该得到一个非常接近0.6617的答案 - 这当然是正确的答案！通过改变传递给的参数samplePoints()，可以解决您喜欢的任何维度的广义问题。