Integracja w Mathematica

Question

Chciałbym uzyskać inne rozwiązanie problemu, który rozwiązałem "symbolicznie" i przez małą symulację. Teraz chciałbym wiedzieć, w jaki sposób mogę uzyskać integrację bezpośrednio za pomocą Mathematica.

Proszę wziąć pod uwagę cel reprezentowany przez dysk z r = 1, wyśrodkowany na (0,0). Chcę zrobić symulację mojego prawdopodobieństwa trafienia w rzutki rzucania celu.

Teraz nie mam żadnych uprzedzeń je rzucają, czyli średnio będę hit centrum il = 0, ale mój wariancja jest 1.

Zważywszy współrzędnych mojego dart, jak trafić w cel (lub ściana :-)) I mają następujące rozkłady 2 Gaussians:

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

z tych 2 rozkładu w środku 0, przy równych wariancji = 1, mój łączny rozkład staje się dwuwymiarowe Gaussa, takie jak:

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

Więc potrzebuję znać moje prawdopodobieństwo, aby trafić cel lub prawdopodobieństwo x^2 + y^2 jest gorsze od 1.

Integracja po transformacji w układzie współrzędnych biegunowych dała mi najpierw moje rozwiązanie: .39. Symulacja potwierdziła go za pomocą:

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

czuję się tam bardziej elegancki sposób, aby rozwiązać ten problem przy użyciu zdolności integracyjnych Mathematica, ale nigdy nie dostał do map eter pracy.

Answer 1

Istnieje kilka sposobów na zrobienie tego.

Najprościej byłoby użyć NIntegrate jak:

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} 

Jest to kolejny sposób na to empirycznie (podobny do powyższego przykładu), ale dużo wolniej niż przy użyciu 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

Wbudowana funkcja NProbability jest również szybki:

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

lub

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

dadzą {0.031, 0.393469}.

Ponieważ suma kwadratów n standardowych normalnych jest rozprowadzany ChiSquare[n], otrzymasz bardziej opływowy wyraz NProbability[z < 1, z \[Distributed] ChiSquareDistribution[2]] gdzie z=x^2+y^2 i x i y są rozproszone NormalDistribution[0,1]. Czas jest taki sam jak powyżej: {0.031, 0.393469}.

EDYCJA: Czasy są dla Vista 64bit Core2 Duo T9600 2.80GHz z 8G pamięci (MMA 8.0.4). Rozwiązanie Yoda na tej maszynie ma czas {0.031, 0.393469}.

EDIT 2: Symulacje pomocą ChiSquareDistribution[2] może być wykonane w następujący sposób:

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

daje {0.031, 0.3946}.

EDIT 3: Więcej szczegółów na temat taktowania:

Dla

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

uzyskać {0.047, 0.031, 0.031, 0.031, 0.031, 0.016, 0.016, 0.031, 0.015, 0.016}

Dla

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

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

Dla

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

uzyskać {0.047, 0.015, 0.016, 0.016, 0.031, 0.015, 0.016, 0.016, 0.015, 0.}.

Dla Yody

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

uzyskać {0.031, 0.032, 0.015, 0., 0.016, 0., 0.015, 0.016, 0.016, 0.}.

Dla empirycznej oceny

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

mam {0.031, 0.016, 0.016, 0., 0.015, 0.016, 0.015, 0., 0.016, 0.016}.