Решите «неразрешимые» с помощью методов Монте-Карло

reshite nerazreshimye s pomoshhyu metodov monte karlo?v=1656636021

Как решить «неразрешимую» проблему?

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

К счастью, существуют методы, которые могут приблизить решение этих проблем с помощью простого приема.

Методы Монте-Карло это класс методов, которые могут быть применены к вычислительно «сложным» проблемам для получения почти достаточно точных ответов. Общая предпосылка очень проста:

  1. Произвольно произведите входные данные для проблемы
  2. Для каждого образца вычислите результат
  3. Объедините результаты, чтобы приблизить решение

Как аналогию, представьте, что вы муравей ползающий по большой плиточной мозаике. С вашей точки зрения, вам нелегко понять, что изображено на мозаике.

Если бы вы начали ходить по мозаике и пробовать посещаемые через случайные промежутки времени плитки, вы бы составили приблизительное представление о том, что мозаика показывает. Чем больше образцов вы возьмете, тем лучше будет ваше приближение.

Если бы вы могли покрыть каждую отдельную плитку, вы бы получили идеальное представление мозаики. Однако в этом нет необходимости – после определенного количества выборок у вас будет достаточно хорошая оценка.

Это как раз то, как методы Монте-Карло приближают решение «неразрешимых» проблем.

Название относится к известному казино в Монако. Его придумал в 1949 году один из пионеров метода Станислав Улам. Как сообщается, дядя Улама был азартным игроком, и связь между элементом «случайности» в азартных играх и методами Монте-Карло, по-видимому, была особенно очевидна Станиславу.

Лучший способ понять техническую концепцию — погрузиться в нее и увидеть, как она работает. Остальная статья покажет, как методы Монте-Карло могут решить три интересных проблемы. Примеры будут представлены на языке программирования Julia.

Знакомство с Юлией

Есть несколько языков, которые вы можете выучить, если вы заинтересованы в специализации в области данных. Одним из вариантов, который в последние годы становится все более серьезным, является речь под названием Джулия.

Julia – это язык численного программирования, принятый в ряде количественных дисциплин. Его можно скачать бесплатно. Существует очень хороший интерфейс на основе браузера под названием JuliaBox, который работает на Jupyter Notebook.

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

Идя параллельно

Чтобы запустить Julia в нескольких процессах, перейдите к терминалу (или откройте новый сеанс терминала в JuliaBox) и выполните следующую команду:

$ julia -p 4

Это инициирует сеанс Julia на четырех ЦБ. Чтобы определить функцию в Julia, используется такой синтаксис:

function square(x) return x^2 end

Правильно — вместо отступлений или фигурных скобок Юлия использует a start-end подход. Циклы for похожие:

for i = 1:10 print(i) end

Конечно, вы можете добавить пробелы и отступы, чтобы облегчить читабельность.

Возможности Джулии относительно параллельного программирования базируются в основном на двух концепциях: удаленные справки и удаленные звонки.

  • Удаленные ссылки – это объекты, которые по существу действуют как именуемые заполнители для объектов, определенных в других процессах.
  • Удаленные вызовы позволяют процессам вызывать функции на основе аргументов, хранящихся в других процессах.

Важно определить функции для всех процессов. Просмотрите код ниже:

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

The @everywhere макрос обеспечивает hello() функция определена для всех процессов. The @spawn макрос используется для обертывания замыкания вокруг выражения hello("World!")затем оценивается удаленно в автоматически выбранном процессе.

Результат этого выражения мгновенно возвращается как a Future дистанционная справка. Если вы попытаетесь напечатать result, вы будете разочарованы. Выход из hello("World!") было оценено в другом процессе, и его здесь нет. Чтобы сделать его доступным, используйте fetch() метод.

Если нерест и выборка кажутся слишком тяжелыми, вам повезло. Юлия также имеет a @parallel макрос, который примет часть тяжелой работы, необходимой для параллельного выполнения задач.

@parallel работает как автономно, так и с функцией «уменьшения» для сбора результатов всех процессов и сведения их к конечному результату. Посмотрите на код ниже:

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

Цикл for просто возвращает значение i на каждом шагу. The @parallel макрос использует оператор добавления в качестве редуктора. Он принимает каждое значение i и прибавляет его к предыдущим значениям.

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

Помня возможности параллельного программирования Джулии, давайте перейдем к тому, как мы можем использовать методы Монте-Карло для решения интересных примеров задач.

Игра в лотерею

В качестве первого примера представим, что мы играем в лотерею. Идея проста – выберите шесть уникальных чисел от 1 до 50. Каждый билет стоит, скажем, 2 фунта.

  • Если вы совпадаете все шесть чисел с вытянутыми, вы выиграете большой приз (£1 000 000)
  • Если вы найдете пять чисел, вы выиграете средний приз (100000 фунтов стерлингов)
  • Если вы совпадаете четыре числа, вы выиграете небольшой приз (£100)
  • Если вы совпадаете три числа, вы выиграете очень маленький приз (£10)

Что бы вы ожидали выиграть, если бы играли в эту лотерею каждый день в течение двадцати лет?

Вы можете решить это с помощью ручки и бумаги, используя немного теории вероятностей. Но это займет много времени! Почему бы не использовать метод Монте-Карло?

Подход почти подозрительно прост – симулируйте игру много раз и усредняйте результат.

Начните Юлия:

$ 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

Количество соответствующих чисел исчисляется с помощью Юлии indexin() функция. Это берет массив и для каждого элемента возвращает индекс его позиции в другом массиве (или ноль, если элемент не найден). В отличие от многих современных языков Julia индексирует от единицы, а не от нуля.

The .!= 0 синтаксис проверяет, какие из этих индексов не равны нулю, и возвращает любой из них true или false для каждого. Наконец, количество true‘s суммируются, давая общие соответствующие числа.

Теперь давайте смоделируем игру в лотерею каждый день в течение двадцати лет… десять тысяч раз параллельно.

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

Небольшая отдача, верно?

Вы можете расширить код, чтобы включить более сложные правила и сценарии, и увидеть, как это повлияет на результат.

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

Давайте посмотрим, что еще позволяют нам сделать методы Монте-Карло.

Значение числа пи

Пи (или π) – это математическая константа. Пожалуй, он наиболее известен своим появлением в формуле для площади круга:

A = πr²

π является примером иррационального числа. Его точное значение невозможно представить в виде какой-либо доли двух целых чисел. На самом деле π также является примером трансцендентного числа — нет даже полиномиальных уравнений, для которых оно является решением.

Вы можете подумать, что это делает получение точного значения для непростым. Или это так?

На самом деле вы можете найти достаточно хорошую оценку π, используя метод Монте-Карло. Визуальная аналогия может выглядеть так:

  • Начертите на стене квадрат 2м×2м. Внутри нарисуйте круг радиусом 1м.
  • Теперь сделайте несколько шагов назад и наугад швырните краской в ​​стену. Считайте каждый раз, когда краска попадает на круг.
  • После ста бросков вычислите, какая часть бросков попала в круг. Умножьте это на площадь квадрата. Вот ваша оценка для π.
Nx1UB51Yol0GpIObWhytD9phm5kkZLEgftGw
Площадь круга составляет примерно 78% площади квадрата, так что примерно такой процент краски попадает внутрь круга

Причина, по которой это работает, очень интуитивно ясна. При выборе случайных точек из квадрата, содержащего круг, вероятность выбора точек из круга пропорциональна площади круга.

Имея достаточное количество случайных выборок, мы можем найти надежную оценку этой доли, с.

Теперь мы знаем, что площадь квадрата равна 2×2 = 4 м², а площадь окружности равна π×р². Поскольку радиус р равна 1, площадь круга равна просто π.

Поскольку мы знаем площадь квадрата и имеем оценку пропорции с его площади, охваченной кругом, можно оценить π. Просто умножьте с×4.

Чем больше случайных выборок мы берем, тем лучше оценка с будет. Однако выигрыш в точности уменьшается, поскольку мы берём все больше и больше образцов.

Вот код Julia для симуляции этого примера. Я запустил это в терминале JuliaBox, используя следующую команду для запуска 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 использует уравнение круга, чтобы проверить, лежат ли точки внутри воображаемой цепи, подсчитывая количество обращений. Функция возвращает долю обращений, умноженную на четыре.

Параллельный запуск этой функции позволит получить очень большое количество образцов, обеспечивая гораздо большую точность.

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

print(Pi)

The nworkers() метод возвращает количество используемых процессоров (в этом случае четыре). Это означает, что каждый процесс выполняется throwPaint() метод, сто миллионов раз. В общем, это дает нам множество образцов и делает очень точную оценку значения π.

Более широкая картина: интеграция

Приведенный выше пример π является конкретным примером более общего использования приближения Монте-Карло — решение проблем интеграции.

Интегрирование – это техника счисления, которая находит область, определенную математической функцией. Например, простая кривая может быть определена функцией:

f(x) → x²

И соответствующий график будет таким:

dD-8jjpkSQZNigDOy-Bn12P1SAD1YGp-0CZn
f(x) → x² дает классическую U-образную кривую, проходящую через начало координат

Площадь под кривой определяется путём интегрирования f(x).

gN77Fve6MF4oTwJhJFPEKXnCh8cPkB0vswMj

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

В малых размерах площадь под кривой может быть аппроксимирована относительно простыми алгоритмами, такими как метод трапеции.

Однако из-за проклятия размерности это становится вычислительно невозможным в высших измерениях. Для оценки площади можно использовать методы на основе Монте-Карло.

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

apgKXkwcEs7p6sDQiFYdth4DAHW2UUh-yNDp
Площадь под кривой составляет около 45% площади квадрата. Вычисления не нужны!

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

Как последний пример возьмем сложную математическую головоломку.

Сложная математическая головоломка

Наугад выберите две точки в единичном кубе. Каково в среднем расстояние между ними?

LycK4mnl7-FALRaqOlQUIYDWAEHcRUTLnV5B
Наугад выберите две точки внутри единичного куба. Каково ожидаемое расстояние между ними?

Сейчас я вас предупрежу математическое решение не совсем тривиальное.

30gI-9USoAyUK1kzjqfli28VwbJ2wpOb04Jl
Если вам нравится решить несколько интегралов, что ж, хорошо для вас! Для остальных из нас всегда Монте-Карло…

Однако можно получить точную оценку с помощью – как вы догадались – метода Монте-Карло.

$ 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

Наконец, запустите эти две функции вместе параллельно. Вместо того чтобы свести к одному результату, на этот раз мы запишем каждый результат в a SharedArray объект. SharedArray объекты позволяют разным процессам получать доступ к данным, хранящимся в одном объекте массива.

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

Вы должны получить ответ, очень близкий к 0,6617 – и это, конечно, правильный ответ! Изменив переданный аргумент samplePoints()Вы можете решить обобщенную задачу в любом количестве измерений.

Что дальше?

Надеемся, вы нашли это вступление к методам Монте-Карло полезным!

Если их правильно внедрить, они предоставляют бесценный инструмент для научных работников по обработке данных, инженеров, финансовых математиков и исследователей… и любого, чья работа предполагает понимание сложных систем.

Если вам интересно узнать больше об их программах, есть масса ресурсов в Интернете. Однако лучший способ обучения – это практика! Когда вы поймете основное положение, почему бы не смоделировать свои примеры Монте-Карло?

Пожалуйста, оставляйте отзывы или комментарии ниже!

Добавить комментарий

Ваш адрес email не будет опубликован.