Integración en Mathematica

Question

Me gustaría obtener una solución diferente a un problema que he resuelto "simbólicamente" y mediante una pequeña simulación. Ahora, me gustaría saber cómo podría haber conseguido la integración directamente usando Mathematica.

Considere un objetivo representado por un disco con r = 1, centrado en (0,0) .Quiero hacer una simulación de mi probabilidad de golpear este objetivo lanzando dardos.

Ahora, no tengo prejuicios tirarlos, es decir, en promedio, voy a golpear el centro de la mu = 0 pero mi varianza es 1.

Teniendo en cuenta las coordenadas de mi dardo al chocar contra el destino (o la pared :-)) tengo las siguientes distribuciones, 2 gaussianas:

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

con los 2 distribución centrada en 0 con igual varianza = 1, mi distribución conjunta se convierte en una gaussiana bivariante como:

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

Así que necesito saber mi probabilidad de alcanzar el objetivo o la probabilidad de que x^2 + y^2 sea inferior a 1.

Una integración después de una transformación en un sistema de coordenadas polares me dio primero mi solución: .39. Simulación confirmó usando:

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

me siento hubo forma más elegante de resolver este problema utilizando las capacidades de integración de Mathematica, pero nunca llegó a mapear el trabajo éter.

Answer 1

En realidad, hay varias formas de hacerlo.

Lo más sencillo sería utilizar NIntegrate como:

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} 

Ésta es otra manera de hacerlo empíricamente (similar a su ejemplo anterior), pero mucho más lento que usar 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

Función incorporada NProbability también es rápido:

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

o

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

ambos dan {0.031, 0.393469}.

Dado que la suma de los cuadrados de n normales estándar se distribuye ChiSquare[n], se obtiene una expresión más simplificada NProbability[z < 1, z \[Distributed] ChiSquareDistribution[2]] donde z=x^2+y^2 y x y y se distribuyen NormalDistribution[0,1]. El tiempo es el mismo que el anterior: {0.031, 0.393469}.

EDITAR: Los tiempos son para una máquina Vista 64bit Core2 Duo T9600 2.80GHz con memoria 8G (MMA 8.0.4). La solución de Yoda en esta máquina tiene sincronización {0.031, 0.393469}.

EDIT 2: Simulación utilizando ChiSquareDistribution[2] se puede hacer como sigue:

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

produce {0.031, 0.3946}.

EDIT 3: Más detalles sobre horarios:

Para

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

me sale {0.047, 0.031, 0.031, 0.031, 0.031, 0.016, 0.016, 0.031, 0.015, 0.016}

Para

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

consigo {0.047, 0.031, 0.032, 0.031, 0.031, 0.016, 0.031, 0.015, 0.016, 0.031}.

Para

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

me sale {0.047, 0.015, 0.016, 0.016, 0.031, 0.015, 0.016, 0.016, 0.015, 0.}.

Para

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

de Yoda consigo {0.031, 0.032, 0.015, 0., 0.016, 0., 0.015, 0.016, 0.016, 0.}.

Para la estimación empírica

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

que tiene {0.031, 0.016, 0.016, 0., 0.015, 0.016, 0.015, 0., 0.016, 0.016}.