6

我想通过一个小模拟得到一个不同的解决方案来解决我“象征性地”解决的问题。现在，我想知道如何直接使用 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 的集成能力有更优雅的方法来解决这个问题，但从来没有映射以太工作。

4

2 回答 2

6

实际上有几种方法可以做到这一点。

最简单的方法是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}
于 2011-12-20T03:38:13.897 回答
4

内置功能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和分布。时间同上：。xyNormalDistribution[0,1]{0.031, 0.393469}

编辑：时序适用于具有 8G 内存 (MMA 8.0.4) 的 Vista 64 位 Core2 Duo T9600 2.80GHz 机器。Yoda 在这台机器上的解决方案有时间{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}。

于 2011-12-20T04:20:59.983 回答