Интеграция в Mathematica

Я хотел бы получить другое решение проблемы, которую я решил «символически» и посредством небольшой симуляции. Теперь я хотел бы знать, как я мог получить интеграцию напрямую с помощью Mathematica.

Пожалуйста, рассмотрите цель, представленную диском с r = 1, с центром в (0,0). Я хочу смоделировать мою вероятность попасть в эту цель, бросая дротики.

Теперь у меня нет смещения, бросающего их, то есть в среднем я попаду в центр mu = 0, но моя дисперсия равна 1.

Учитывая координату моего дротика при попадании в цель (или в стену :-)), у меня есть следующие распределения, 2 гауссиана:

XDistribution : 1/Sqrt[2 \[Pi]\[Sigma]^2] E^(-x^2/(2 \[Sigma]^2))

YDistribution : 1/Sqrt[2 \[Pi]\[Sigma]^2] E^(-y^2/(2 \[Sigma]^2))

С этими 2 распределениями с центром в 0 с равной дисперсией =1 мое совместное распределение становится двумерным гауссовым, например:

1/(2 \[Pi]\[Sigma]^2) E^(-((x^2 + y^2)/(2 \[Sigma]^2)))

Поэтому мне нужно знать мою вероятность попасть в цель или вероятность того, что x ^ 2 + y ^ 2 будет меньше 1.

Интегрирование после преобразования в полярной системе координат дало мне первое решение: .39. Моделирование подтвердило это с помощью:

Total@ParallelTable[
   If[
      EuclideanDistance[{
                         RandomVariate[NormalDistribution[0, Sqrt[1]]], 
                         RandomVariate[NormalDistribution[0, Sqrt[1]]]
                        }, {0, 0}] < 1, 1,0], {1000000}]/1000000

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


integration gaussian wolfram-mathematica
person 500    schedule 20.12.2011    source источник

Ответы (2)


arrow_upward
6
arrow_downward

На самом деле есть несколько способов сделать это.

Проще всего было бы использовать NIntegrate как:

JointDistrbution = 1/(2 \[Pi] \[Sigma]^2) E^(-((x^2 + y^2)/(2 \[Sigma]^2)));
NIntegrate[JointDistrbution /. \[Sigma] -> 1, {y, -1, 1}, 
    {x, -Sqrt[1 - y^2], Sqrt[1 - y^2]}] // Timing

Out[1]= {0.009625, 0.393469}

Это еще один способ сделать это эмпирически (аналогично вашему примеру выше), но намного медленнее, чем использование NIntegrate:

(EuclideanDistance[#, {0, 0}] <= 1 & /@ # // Boole // Total)/
     Length@# &@RandomVariate[NormalDistribution[0, 1], {10^6, 2}] // 
  N // Timing

Out[2]= {5.03216, 0.39281}
person abcd    schedule 20.12.2011
comment
Мне показалось интересным, что Mathematica также может использовать Integrate[] JointDistribution. - person Daniel Chisholm; 21.12.2011

arrow_upward
4
arrow_downward

Встроенная функция NProbability тоже быстрая:

NProbability[ x^2 + y^2 <= 1, {x, y} \[Distributed] 
BinormalDistribution[{0, 0}, {1, 1}, 0]] // Timing

or

NProbability[x^2 + y^2 <= 1, x \[Distributed] 
NormalDistribution[0, 1] && y \[Distributed] 
NormalDistribution[0, 1] ] // Timing

оба дают {0.031, 0.393469}.

Поскольку сумма квадратов n стандартных нормалей распределяется ChiSquare[n], вы получаете более простое выражение NProbability[z < 1, z \[Distributed] ChiSquareDistribution[2]], где z=x^2+y^2, x и y распределяются NormalDistribution[0,1]. Время такое же, как указано выше: {0.031, 0.393469}.

РЕДАКТИРОВАТЬ: Время указано для машины Vista 64bit Core2 Duo T9600 2,80 ГГц с памятью 8G (MMA 8.0.4). Решение Йоды на этой машине имеет время {0.031, 0.393469}.

РЕДАКТИРОВАТЬ 2: Моделирование с использованием ChiSquareDistribution[2] можно выполнить следующим образом:

(data = RandomVariate[ChiSquareDistribution[2], 10^5]; 
  Probability[w <= 1, w \[Distributed] data] // N) // Timing

дает {0.031, 0.3946}.

РЕДАКТИРОВАТЬ 3: Подробнее о таймингах:

За

First@Transpose@Table[Timing@
  NProbability[x^2 + y^2 <= 1, {x, y} \[Distributed] 
  BinormalDistribution[{0, 0}, {1, 1}, 0]], {10}]

я получаю {0.047, 0.031, 0.031, 0.031, 0.031, 0.016, 0.016, 0.031, 0.015, 0.016}

За

First@Transpose@Table[Timing@
NProbability[x^2 + y^2 <= 1, 
 x \[Distributed] NormalDistribution[0, 1] && 
  y \[Distributed] NormalDistribution[0, 1] ], {10}]

Я получаю {0.047, 0.031, 0.032, 0.031, 0.031, 0.016, 0.031, 0.015, 0.016, 0.031}.

За

First@Transpose@Table[Timing@
NProbability[z < 1, 
 z \[Distributed] ChiSquareDistribution[2]], {10}]

Я получаю {0.047, 0.015, 0.016, 0.016, 0.031, 0.015, 0.016, 0.016, 0.015, 0.}.

Для йоды

First@Transpose@Table[Timing@(JointDistrbution = 
  1/(2 \[Pi] \[Sigma]^2) E^(-((x^2 + y^2)/(2 \[Sigma]^2))); 
 NIntegrate[
  JointDistrbution /. \[Sigma] -> 1, {y, -1, 
   1}, {x, -Sqrt[1 - y^2], Sqrt[1 - y^2]}]), {10}]

Я получаю {0.031, 0.032, 0.015, 0., 0.016, 0., 0.015, 0.016, 0.016, 0.}.

Для эмпирической оценки

First@Transpose@Table[Timing@(Probability[w <= 1, 
 w \[Distributed] data] // N), {10}]

У меня {0.031, 0.016, 0.016, 0., 0.015, 0.016, 0.015, 0., 0.016, 0.016}.

person kglr    schedule 20.12.2011
comment
Я нахожу очень подозрительным, что ваши тайминги точно одинаковы для всех трех ваших решений и моего... Я, конечно, получаю очень разные тайминги - person abcd; 20.12.2011
comment
@yoda, любопытно, не так ли? Я собирался спросить вас, можете ли вы запустить приведенный выше код на своем компьютере. - person kglr; 20.12.2011
comment
Это время, которое я получаю для каждого из трех ваших методов (в порядке, который вы указали) и моего (последнего): {0.035673, 0.022273, 0.097494, 0.009067} - person abcd; 20.12.2011
comment
@йода, спасибо. Я вижу много вариаций (и несколько повторяющихся чисел) в своих таймингах. Очевидно, что-то неслучайное в моем окружении? - person kglr; 20.12.2011
comment
Кстати, я попробовал те же самые тайминги со свежим сеансом ядра и 100 повторениями, и я получил аналогичный шаблон циклирования. - person kglr; 20.12.2011
comment
Хм... понятия не имею, почему! Спасибо хоть :) - person abcd; 20.12.2011