Scala中Pi的蒙特卡罗计算答案

【问题标题】：Monte Carlo calculation of Pi in ScalaScala中Pi的蒙特卡罗计算
【发布时间】：2014-09-04 21:17:58
【问题描述】：

假设我想用蒙特卡罗模拟计算 Pi 作为练习。

我正在编写一个函数，它在正方形(0, 1), (1, 0) 中随机选取一个点并测试该点是否在圆内。

import scala.math._
import scala.util.Random

def circleTest() = {
  val (x, y) = (Random.nextDouble, Random.nextDouble)
  sqrt(x*x + y*y) <= 1
}

然后我正在编写一个函数，它将测试函数和试验次数作为参数，并返回发现测试为真的试验分数。

def monteCarlo(trials: Int, test: () => Boolean) =
  (1 to trials).map(_ => if (test()) 1 else 0).sum * 1.0 / trials

...我可以计算 Pi

monteCarlo(100000, circleTest) * 4

现在我想知道monteCarlo功能是否可以改进。你如何写monteCarlo 高效和可读？

例如，由于试验次数很多，是否值得使用 view 或 iterator 代替 Range(1, trials) 和 reduce 代替 map 和 sum ？

【问题讨论】：

主要的加速应该来自删除 sqrt。 sqrt(x*x + y*y) <= 1 与 x*x + y*y <= 1 相同（即两边平方）
好点！谢谢。
我在这里放了一个包含各种方法的微基准测试页面：gist.github.com/willf/e1f2ce95f04442af53e5 我很高兴地说递归是最快的:)
@WillFitzgerald，很有趣，谢谢。我实际上很惊讶基于 Stream 的版本有多快。在这种事情上，递归可能总是最快的，因为开销较小，但是，正如他们所说，“递归是函数式编程的首选”——高效，但对于工作程序来说不一定是最清晰或最快的。跨度>
在优化时，想想它编译成机器码的样子。机器代码中最快的事情总是简单的 while 循环。方法调用、对象分配（对于 lambdas 通常是必需的）等在像这样的紧密循环中非常昂贵。所以最快的将是一个while循环。下面的递归方法只是因为尾调用优化而很快，这意味着它们不会编译为字节码中的递归方法调用，而是编译为 while 循环。

标签： scala montecarlo

【解决方案1】：

值得注意的是，Random.nextDouble 是有副作用的——当你调用它时，它会改变随机数生成器的状态。这可能不是您关心的问题，但由于这里已经有五个答案，我认为添加一个纯功能性的答案不会有什么坏处。

首先，您需要一个随机数生成单子实现。幸运的是，NICTA 提供了与 Scalaz 集成的a really nice one。你可以这样使用它：

import com.nicta.rng._, scalaz._, Scalaz._

val pointInUnitSquare = Rng.choosedouble(0.0, 1.0) zip Rng.choosedouble(0.0, 1.0)

val insideCircle = pointInUnitSquare.map { case (x, y) => x * x + y * y <= 1 }

def mcPi(trials: Int): Rng[Double] =
  EphemeralStream.range(0, trials).foldLeftM(0) {
    case (acc, _) => insideCircle.map(_.fold(1, 0) + acc)
  }.map(_ / trials.toDouble * 4)

然后：

scala> val choosePi = mcPi(10000000)
choosePi: com.nicta.rng.Rng[Double] = com.nicta.rng.Rng$$anon$3@16dd554f

尚未计算任何内容——我们刚刚建立了一个计算，该计算将在执行时随机生成我们的值。为方便起见，让我们在IO monad 中就地执行它：

scala> choosePi.run.unsafePerformIO
res0: Double = 3.1415628

这不会是最高效的解决方案，但它已经足够好了，对于许多应用程序来说可能不是问题，并且引用透明性可能是值得的。

【讨论】：

非常感谢！我在考虑Random 的副作用以及如何使蒙特卡洛变得纯净。我不接受这个答案，因为它超出了问题范围，但我很感激并开始学习 Random monad。

【解决方案2】：

基于流的版本，另一种选择。我认为这很清楚。

def monteCarlo(trials: Int, test: () => Boolean) =
    Stream
      .continually(if (test()) 1.0 else 0.0)
      .take(trials)
      .sum / trials

（sum 不是专门用于流的，但实现（在 TraversableOnce 中）只是调用专门的 foldLeft 并且“允许 GC 沿途收集。”因此 .sum 不会强制流进行评估，因此不会一次将所有试验都保存在内存中）

【讨论】：

有趣的是，reduceLeft 以这种方式专业化。我喜欢这个解决方案，但我认为reduceLeft 是专业的这一事实使它不那么直观。
不确定你的意思。 reduceLeft 在任何地方都可以工作并且做同样的事情，只是流版本有一个特定的效率实现。许多其他集合方法也是如此——它们可能对某些集合类型有专门的实现。你可以用 sum 替换这个 reduceLeft，它只会强制创建流内容并将其保存在内存中（我认为，没有检查）
啊，不，只是查看了源代码，而 .sum 只是调用了 foldLeft，所以在流上也应该很有效。我会编辑我的答案
当然，我不想将流的所有 10000 个元素都保存在内存中。所以从性能的角度来看，reduceLeft 很好。另一方面，恐怕人们只是不知道 reduceLeft 的这种特殊行为。所以他们会认为流在内存中保存了它的所有元素。
@WillFitzgerald 感谢您抽出时间对其进行基准测试！

【解决方案3】：

我认为以下递归版本没有问题：

def monteCarlo(trials: Int, test: () => Boolean) = {
  def bool2double(b: Boolean) = if (b) 1.0d else 0.0d
  @scala.annotation.tailrec
  def recurse(n: Int, sum: Double): Double = 
    if (n <= 0) sum / trials
    else recurse(n - 1, sum + bool2double(test()))
  recurse(trials, 0.0d)
}

【讨论】：

【解决方案4】：

还有一个 foldLeft 版本：

def monteCarloFold(trials: Int, test: () => Boolean) = 
  (1 to trials).foldLeft(0.0d)((s,i) => s + (if (test()) 1.0d else 0.0d)) / trials

这比问题中的map 版本更节省内存。

【讨论】：

【解决方案5】：

使用尾递归可能是一个想法：

def recMonteCarlo(trials: Int, currentSum: Double, test:() => Boolean):Double = trials match {
  case 0 => currentSum
  case x => 
    val nextSum = currentSum + (if (test()) 1.0 else 0.0)
    recMonteCarlo(trials-1, nextSum, test)

def monteCarlo(trials: Int, test:() => Boolean) = {
  val monteSum = recMonteCarlo(trials, 0, test)
  monteSum / trials
}

【讨论】：

谢谢。虽然它有点冗长，不是吗？
优点是您不必为了计算总和而将整个列表保存在内存中。这意味着您可以更轻松地扩展此版本。
我明白了。不过，我希望还有其他解决方案，它们也可以节省内存。
当然。总是有不止一种方法:)

【解决方案6】：

在并行集合上使用aggregate，像这样，

def monteCarlo(trials: Int, test: () => Boolean) = {
  val pr = (1 to trials).par
  val s = pr.aggregate(0)( (a,_) => a + (if (test()) 1 else 0), _ + _) 
  s * 4.0 / trials
}

部分结果与其他测试计算并行汇总。

【讨论】：

谢谢。并行性很有趣，但我现在不想处理它。