# 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]:= XPolygonRegionFunction[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]:= XPolygonProbability[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]:= XRegularPolygonProbability[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]:= XWith[{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]=