New in Wolfram
Mathematica
8: New and Improved Core Algorithms
◄
previous

next
►
Core Algorithms
Compute Binormal Probability over a Polygon
Probability of standard binormal variate with
over a regular octagon of unit radius expressed in terms of Owen's T function.
In[1]:=
X
PolygonRegionFunction[pts_] := Block[{ppts = Join[pts, pts], diff, x, y}, Function[{x, y}, Evaluate[ And @@ Table[diff = ppts[[i + 1]]  ppts[[i]]; Cross[diff].({x, y}  ppts[[i]]) >= 0, {i, Length[pts]}]]]]
In[2]:=
X
PolygonProbability[pts : {{_, _} ..}, BinormalDistribution[{mu1_, mu2_}, {s1_, s2_}, r_]? DistributionParameterQ] /; PolygonRegionFunction[pts][mu1, mu2] := Module[{tpts = Join[pts, pts[[{1}]]], res = 0, h1, k1, h2, k2, cross, diff, h, a1, a2}, tpts = Transpose[{{1, 1}/Sqrt[2 (1 + r)], {1, 1}/Sqrt[ 2 (1  r)]}.((Transpose[pts]  {mu1, mu2})/{s1, s2})]; Do[{{h1, k1}, {h2, k2}} = tpts[[{i, i + 1}]]; cross = Abs[h1 k2  k1 h2]; If[! PossibleZeroQ[cross], diff = {h2  h1, k2  k1}; h = cross/Sqrt[diff.diff]; a1 = {h1, k1}.diff/cross; a2 = {h2, k2}.diff/cross; res += OwenT[h, Max[a1, a2]]  OwenT[h, Min[a1, a2]]], {i, Length[pts]}]; res]
In[3]:=
X
RegularPolygonProbability[sides_Integer /; sides >= 3, r_] := 1  FullSimplify[ PolygonProbability[ Most[Table[{Cos[x], Sin[x]}, {x, 0, 2 \[Pi], (2 \[Pi])/sides}]], BinormalDistribution[{0, 0}, {1, 1}, r]], 1 < r < 1]
In[4]:=
X
With[{sides = 8, r = 1/3}, Show[Plot3D[ PDF[BinormalDistribution[r], {x, y}], {x, 2, 2}, {y, 2, 2}, PlotStyle > Opacity[0.5, Hue[.15, 0.2, .8]], Filling > Axis, FillingStyle > {Black, Orange}, Lighting > "Neutral", BoundaryStyle > None, RegionFunction > PolygonRegionFunction[ Most[Table[{Cos[x], Sin[x]}, {x, 0, 2 \[Pi], (2 \[Pi])/sides}]]], Mesh > None, PlotRange > All, NormalsFunction > None, PlotPoints > 50], Plot3D[PDF[BinormalDistribution[r], {x, y}], {x, 3, 3}, {y, 3, 3}, PlotStyle > Opacity[0.2], Mesh > 6, MeshStyle > Gray, PlotRange > All, PlotPoints > 50], Boxed > False, ImageSize > 580, Epilog > Inset[Framed[ Pane[Style[ Row[{TraditionalForm[Pr[1/3]] " \[Equal] ", 1, FullSimplify[4 OwenT[Sqrt[1/2 + 1/Sqrt[2]]/Sqrt[ Sqrt[2]  r], (Sqrt[2]  1) Sqrt[(1 + r)/(1  r)]]  4 OwenT[Sqrt[1/2 + 1/Sqrt[2]]/Sqrt[Sqrt[2]  r], ( Sqrt[2]  1  r)/Sqrt[1  r^2]]  4 OwenT[Sqrt[1/2 + 1/Sqrt[2]]/Sqrt[ Sqrt[2] + r], (Sqrt[2]  1) Sqrt[(1  r)/(1 + r)]]  4 OwenT[Sqrt[1/2 + 1/Sqrt[2]]/Sqrt[Sqrt[2] + r], ( Sqrt[2]  1 + r)/Sqrt[1  r^2]]]}], 11.5], {500, 118}], Background > LightYellow], {Center, Top}, {Center, Top}]]]
Out[4]=