\documentclass{article}

\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{nicefrac}

\begin{document}
\section*{Bayesian Regression Model}

\section*{Beta Mode-Concentration Bayesian Estimator}

We will apply polar coordinates $(r,t)$ on a shifted origin:
$\mathbf{1}=\left[\begin{matrix}1\\1\end{matrix}\right]$.
First, all random variables in the formula will be translated by the origin
$\mathbf{1}$, but importantly, the origin does not need to match
$E[\mathbf{v}]$ (the likelihood is not centered).
The rays extending from the origin $\mathbf{1}$ each have a unique value for the Beta
distribution mode (denoted $\omega$) in the positive quarter-plane where the
mode is defined.
Integrating $r$ and $t$ provides one method to calculate the following normal
CDF for the positive quarter-plane region: $P(\alpha\ge 1, \beta\ge 1)$.
However, sophisticated numerical integration methods already exist.
The benefit of this method is that for any origin, the inner integral for one
ray has an analytical solution.
Therefore, for the custom inference here, only a single (rather than double)
numerical integration is required.

\begin{equation*}
  \begin{aligned}
    \mathbf{v} &= \theta \left[ \begin{matrix} \sin t \\ \cos t \end{matrix} \right] \\
    \omega &= \frac{\sin t}{\sin t+\cos t} \\
  \end{aligned}
\end{equation*}
\begin{equation*}
  \begin{aligned}
    & E\left[ \omega | \alpha\ge 1, \beta\ge 1 \right] \cdot P(\alpha\ge 1, \beta\ge 1) \\
    % Line 2: vector v
    &=
    \frac{1}{2\pi \sqrt{\text{det }\Sigma}}
    \int_0^{\nicefrac{\pi}{2}}
    \omega
    \int_0^\infty
    e^{-\frac{1}{2}(\mathbf{v}-(\mathbf{\mu}-\mathbf{1}))^\intercal \Sigma^{-1} (\mathbf{v}-(\mathbf{\mu}-\mathbf{1}))}
    \theta
    \text{d}\theta
    \text{d}t
    \\
    &=
    % Line 3: Expand the subtraction in the quadratic form.
    \frac{1}{2\pi \sqrt{\text{det }\Sigma}}
    \int_0^{\nicefrac{\pi}{2}}
    \omega
    \int_0^\infty
    e^{
      -\frac{1}{2} \mathbf{v}^\intercal \Sigma^{-1} \mathbf{v}
      + \mathbf{v}^\intercal \Sigma^{-1} (\mathbf{\mu}-\mathbf{1})
      -\frac{1}{2}(\mathbf{\mu}-\mathbf{1})^\intercal \Sigma^{-1} (\mathbf{\mu}-\mathbf{1})
    }
    \theta
    \text{d}\theta
    \text{d}t
    \\
    % Line 4: Substitute theta [sin t cos t] for vector v
    &=
    \frac{1}{2\pi \sqrt{\text{det }\Sigma}}
    \int_0^{\nicefrac{\pi}{2}}
    \omega
    \int_0^\infty
    e^{
      -\frac{1}{2}\theta^2 \left[\begin{matrix} \sin t \\ \cos t \end{matrix}\right]^\intercal
      \Sigma^{-1}
      \left[\begin{matrix} \sin t \\ \cos t \end{matrix}\right]
      + \theta \left[\begin{matrix} \sin t \\ \cos t \end{matrix}\right]^\intercal
      \Sigma^{-1} (\mathbf{\mu}-\mathbf{1})
      -\frac{1}{2} (\mathbf{\mu}-\mathbf{1})^\intercal \Sigma^{-1} (\mathbf{\mu}-\mathbf{1})
    }
    \theta
    \text{d}\theta
    \text{d}t
    \\
    % Line 5: Substitute a, b, c (coefficients that do not depend on theta)
    &=
    \frac{1}{2\pi \sqrt{\text{det }\Sigma}}
    \int_0^{\nicefrac{\pi}{2}}
    \omega
    \int_0^\infty
    e^{
      -\frac{1}{2}\left(a \theta^2
      - 2b \theta
      +c\right)
    }
    \theta
    \text{d}\theta
    \text{d}t
    \\
    % Line 6: Complete the square with a and b.
    &=
    \frac{1}{2\pi \sqrt{\text{det }\Sigma}}
    \int_0^{\nicefrac{\pi}{2}}
    \omega
    e^{\frac{b^2}{2a}-\frac{1}{2}c}
    \int_0^\infty
    e^{
      -\frac{a}{2}\left(\theta^2
      - 2\frac{b}{a} \theta
      + \frac{b^2}{a^2}\right)
    }
    \theta
    \text{d}\theta
    \text{d}t
    \\
    % Line 7: Get ready for u-substitution (translating variable theta by a constant)
    &=
    \frac{1}{2\pi \sqrt{\text{det }\Sigma}}
    \int_0^{\nicefrac{\pi}{2}}
    \omega
    e^{\frac{b^2}{2a}-\frac{1}{2}c}
    \int_0^\infty
    e^{-\frac{a}{2}(\theta-\nicefrac{b}{a})^2}
    \left(\theta - \nicefrac{b}{a}\right)
    +
    \frac{b}{a} e^{-\frac{a}{2}(\theta-\nicefrac{b}{a})^2}
    \text{d}\theta
    \text{d}t
    \\
    % Line 8: Apply the known antiderivatives for x exp(-x^2/2) and exp(-x^2/2).
    &=
    \frac{1}{2\pi \sqrt{\text{det }\Sigma}}
    \int_0^{\nicefrac{\pi}{2}}
    \omega
    e^{\frac{b^2}{2a}-\frac{1}{2}c}
    \left(
      \frac{1}{a}
      e^{-\nicefrac{b^2}{2a^3}}
      +
      \frac{b}{a}
      \sqrt{
      \frac{2\pi}{a}
      }
      F_{\mathcal{N}}\left(\nicefrac{b}{\sqrt{a}}\right)
    \right)
    \text{d}t
  \end{aligned}
\end{equation*}

This is an efficient numerical method for calculating the expectation (with
normal likelihood) of an arbitrary random variable $\omega$, which, conditioned
on the angle formed by $\mathbf{v}$ and an arbitrary origin (e.g. $\mathbf{1}$),
is constant. If not an integral from $0$ to $2\pi$, then the result needs to be
divided by the total probability in the region being integrated.

\end{document}