Integrazione in Mathematica

Question

Mi piacerebbe avere una soluzione diversa a un problema che ho risolto "simbolicamente" e attraverso una piccola simulazione. Ora, vorrei sapere come ho potuto ottenere l'integrazione direttamente con Mathematica.

Si prega di considerare un bersaglio rappresentato da un disco con r = 1, centrato a (0,0). Voglio fare una simulazione della mia probabilità di colpire questo bersaglio lanciando frecce.

Ora, non ho pregiudizi loro lancio, che è in media sarò colpire il centro mu = 0 ma il mio varianza è 1.

Considerando le coordinate del mio dardo come colpire il bersaglio (o la parete :-)) ho le seguenti distribuzioni, 2 gaussiane:

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 quei 2 distribuzione centrata a 0 con uguale varianza = 1, la distribuzione congiunta diventa una gaussiana bivariato quali:

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

Quindi ho bisogno di sapere la mia probabilità di colpire il bersaglio o la probabilità di x^2 + y^2 di essere inferiore a 1.

Un'integrazione dopo una trasformazione in un sistema di coordinate polari mi ha dato prima la mia soluzione: .39. Simulazione confermato utilizzando:

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

mi sento c'erano modo più elegante per risolvere questo problema utilizzando le capacità di integrazione di Mathematica, ma non ha mai avuto per mappare il lavoro etere.

Answer 1

In realtà ci sono diversi modi per farlo.

La più semplice sarebbe quella di utilizzare NIntegrate come:

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} 

Questo è un altro modo per farlo empiricamente (simile al vostro esempio di cui sopra), ma molto più lento rispetto all'utilizzo di 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

La funzione integrata NProbability è anche veloce:

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 

sia dare {0.031, 0.393469}.

Dal momento che la somma dei quadrati delle n normali standard, è distribuito ChiSquare[n], si ottiene un'espressione più snella NProbability[z < 1, z \[Distributed] ChiSquareDistribution[2]] dove z=x^2+y^2 e x e y sono distribuiti NormalDistribution[0,1]. La tempistica è la stessa di sopra: {0.031, 0.393469}.

MODIFICA: I tempi sono per una macchina Vista 64bit Core2 Duo T9600 a 2,80 GHz con memoria 8G (MMA 8.0.4). La soluzione di Yoda su questa macchina ha i tempi {0.031, 0.393469}.

EDIT 2: simulazione utilizzando ChiSquareDistribution[2] può essere fatto come segue:

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

cede {0.031, 0.3946}.

EDIT 3: Maggiori dettagli sul timing:

Per

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

ottengo {0.047, 0.031, 0.031, 0.031, 0.031, 0.016, 0.016, 0.031, 0.015, 0.016}

Per

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

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

Per

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

ottengo {0.047, 0.015, 0.016, 0.016, 0.031, 0.015, 0.016, 0.016, 0.015, 0.}.

Per di Yoda

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

ottengo {0.031, 0.032, 0.015, 0., 0.016, 0., 0.015, 0.016, 0.016, 0.}.

Per la stima empirica

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

ho ottenuto {0.031, 0.016, 0.016, 0., 0.015, 0.016, 0.015, 0., 0.016, 0.016}.