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

Question

Я хотел бы получить другое решение проблемы, которую я решил "символически" и с помощью небольшой симуляции.Теперь я хотел бы знать, как я мог получить интеграцию напрямую с помощью 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))

С этими двумя распределениями, центрированными в 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, но так и не удалось отобразить работу эфира.

Answer 1

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

Самым простым было бы использовать 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}

Answer 2

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

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

или

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 64-битный 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}.