Монте-Карло расчет Пи в Скала

Question

Предположим, я хотел бы рассчитать Пи с помощью симуляции Монте-Карло в качестве упражнения.

Я пишу функцию, которая выбирает точку в квадрате (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

... и я могу рассчитать Пи

monteCarlo(100000, circleTest) * 4

Теперь мне интересно, если monteCarlo функция может быть улучшена. Как бы вы написали monteCarlo эффективный и читаемый?

Например, поскольку количество испытаний велико, стоит ли view или же iterator вместо Range(1, trials) а также reduce вместо map а также sum?

Answer 1

Потоковая версия, для другой альтернативы. Я думаю, что это совершенно ясно.

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

(sum не специализируется на потоках, но реализация (в TraversableOnce) просто вызывает foldLeft это специализировано и "позволяет GC собирать по пути". Таким образом,.sum не заставит поток быть оцененным и, следовательно, не сохранит все испытания в памяти сразу)

Answer 2

Стоит отметить, что Random.nextDouble побочный эффект - когда вы вызываете его, он изменяет состояние генератора случайных чисел. Это может не беспокоить вас, но, поскольку здесь уже есть пять ответов, я полагаю, что ничего не помешает добавить чисто функциональный ответ.

Сначала вам понадобится реализация монады генерации случайных чисел. К счастью, NICTA предлагает действительно хороший вариант, интегрированный со Scalaz. Вы можете использовать это так:

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 Монада ради удобства:

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

Это будет не самое эффективное решение, но оно достаточно хорошее, чтобы не создавать проблем для многих приложений, и ссылочная прозрачность может стоить того.

Answer 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)
}

Answer 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 версия в вопросе.

Answer 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
}

Answer 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
}

где частичные результаты суммируются параллельно с другими тестовыми расчетами.